1
/* mpf_div -- Divide two floats.
3
Copyright 1993, 1994, 1996, 2000, 2001 Free Software Foundation, Inc.
5
This file is part of the GNU MP Library.
7
The GNU MP Library is free software; you can redistribute it and/or modify
8
it under the terms of the GNU Lesser General Public License as published by
9
the Free Software Foundation; either version 2.1 of the License, or (at your
10
option) any later version.
12
The GNU MP Library is distributed in the hope that it will be useful, but
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15
License for more details.
17
You should have received a copy of the GNU Lesser General Public License
18
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20
MA 02111-1307, USA. */
27
mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
31
mp_size_t usize, vsize;
32
mp_size_t rsize, tsize;
33
mp_size_t sign_quotient;
41
sign_quotient = usize ^ vsize;
57
rexp = u->_mp_exp - v->_mp_exp;
70
tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB);
80
MPN_ZERO (tp, tsize - usize);
81
rtp = tp + (tsize - usize);
85
/* Normalize the divisor and the dividend. */
86
if (! (vp[vsize-1] & MP_LIMB_T_HIGHBIT))
90
unsigned normalization_steps;
92
count_leading_zeros (normalization_steps, vp[vsize - 1]);
94
/* Shift up the divisor setting the most significant bit of
95
the most significant limb. Use temporary storage not to clobber
96
the original contents of the divisor. */
97
tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
98
mpn_lshift (tmp, vp, vsize, normalization_steps);
101
/* Shift up the dividend, possibly introducing a new most
102
significant word. Move the shifted dividend in the remainder
104
nlimb = mpn_lshift (rtp, up, usize, normalization_steps);
114
/* The divisor is already normalized, as required.
115
Copy it to temporary space if it overlaps with the quotient. */
116
if (vp - rp <= tsize - vsize)
118
mp_ptr tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
119
MPN_COPY (tmp, vp, vsize);
120
vp = (mp_srcptr) tmp;
123
/* Move the dividend to the remainder. */
124
MPN_COPY (rtp, up, usize);
127
q_limb = mpn_divmod (rp, tp, tsize, vp, vsize);
128
rsize = tsize - vsize;
136
r->_mp_size = sign_quotient >= 0 ? rsize : -rsize;