17
18
You should have received a copy of the GNU Lesser General Public License
18
19
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. */
20
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21
MA 02110-1301, USA. */
23
#include <stdio.h> /* for NULL */
23
25
#include "gmp-impl.h"
24
26
#include "longlong.h"
27
1. Develop >= n bits of src.num / src.den, where n is the number of bits
28
in a double. This (partial) division will use all bits from the
30
2. Use the remainder to determine how to round the result.
31
3. Assign the integral result to a temporary double.
32
4. Scale the temporary double, and return the result.
34
An alternative algorithm, that would be faster:
29
/* All that's needed is to get the high 53 bits of the quotient num/den,
30
rounded towards zero. More than 53 bits is fine, any excess is ignored
33
N_QLIMBS is how many quotient limbs we need to satisfy the mantissa of a
34
double, assuming the highest of those limbs is non-zero. The target
35
qsize for mpn_tdiv_qr is then 1 more than this, since that function may
36
give a zero in the high limb (and non-zero in the second highest).
38
The use of 8*sizeof(double) in N_QLIMBS is an overestimate of the
39
mantissa bits, but it gets the same result as the true value (53 or 48 or
40
whatever) when rounded up to a multiple of GMP_NUMB_BITS, for non-nails.
44
Use the true mantissa size in the N_QLIMBS formala, to save a divide step
47
Examine the high limbs of num and den to see if the highest 1 bit of the
48
quotient will fall high enough that just N_QLIMBS-1 limbs is enough to
49
get the necessary bits, thereby saving a division step.
51
Bit shift either num or den to arrange for the above condition on the
52
high 1 bit of the quotient, to save a division step always. A shift to
53
save a division step is definitely worthwhile with mpn_tdiv_qr, though we
54
may want to reassess this on big num/den when a quotient-only division
57
Maybe we could estimate the final exponent using nsize-dsize (and
58
possibly the high limbs of num and den), so as to detect overflow and
59
return infinity or zero quickly. Overflow is never very helpful to an
60
application, and can therefore probably be regarded as abnormal, but we
61
may still like to optimize it if the conditions are easy. (This would
62
only be for float formats we know, unknown formats are not important and
63
can be left to mpn_get_d.)
67
If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
68
padding n with zeros in temporary space.
70
If/when a quotient-only division exists it can be used here immediately.
71
remp is only to satisfy mpn_tdiv_qr, the remainder is not used.
75
An alternative algorithm, that may be faster:
35
76
0. Let n be somewhat larger than the number of significant bits in a double.
36
77
1. Extract the most significant n bits of the denominator, and an equal
37
78
number of bits from the numerator.
42
83
we are done. If they are different, repeat the algorithm from step 1,
43
84
but first let n = n * 2.
44
85
4. If we end up using all bits from the numerator and denominator, fall
45
back to the first algorithm above.
86
back to a plain division.
46
87
5. Just to make life harder, The computation of a + 1 and b + 1 above
47
88
might give carry-out... Needs special handling. It might work to
48
89
subtract 1 in both cases instead.
91
Not certain if this approach would be faster than a quotient-only
92
division. Presumably such optimizations are the sort of thing we would
93
like to have helping everywhere that uses a quotient-only division. */
52
96
mpq_get_d (const MP_RAT *src)
56
101
mp_size_t nsize = src->_mp_num._mp_size;
57
102
mp_size_t dsize = src->_mp_den._mp_size;
58
mp_size_t qsize, rsize;
59
mp_size_t sign_quotient = nsize ^ dsize;
103
mp_size_t qsize, prospective_qsize, zeros, chop, tsize;
104
mp_size_t sign_quotient = nsize;
61
106
#define N_QLIMBS (1 + (sizeof (double) + BYTES_PER_MP_LIMB-1) / BYTES_PER_MP_LIMB)
62
107
mp_limb_t qarr[N_QLIMBS + 1];
111
ASSERT (dsize > 0); /* canonical src */
113
/* mpn_get_d below requires a non-zero operand */
114
if (UNLIKELY (nsize == 0))
70
118
nsize = ABS (nsize);
71
119
dsize = ABS (dsize);
72
120
np = src->_mp_num._mp_d;
73
121
dp = src->_mp_den._mp_d;
75
rsize = dsize + N_QLIMBS;
76
rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB);
78
/* Normalize the denominator, i.e. make its most significant bit set by
79
shifting it NORMALIZATION_STEPS bits to the left. Also shift the
80
numerator the same number of steps (to keep the quotient the same!). */
81
if ((dp[dsize - 1] & GMP_NUMB_HIGHBIT) == 0)
123
prospective_qsize = nsize - dsize + 1; /* from using given n,d */
124
qsize = N_QLIMBS + 1; /* desired qsize */
126
zeros = qsize - prospective_qsize; /* padding n to get qsize */
127
exp = (long) -zeros * GMP_NUMB_BITS; /* relative to low of qp */
129
chop = MAX (-zeros, 0); /* negative zeros means shorten n */
132
zeros += chop; /* now zeros >= 0 */
134
tsize = nsize + zeros; /* size for possible copy of n */
85
unsigned normalization_steps;
87
count_leading_zeros (normalization_steps, dp[dsize - 1]);
88
normalization_steps -= GMP_NAIL_BITS;
90
/* Shift up the denominator setting the most significant bit of
91
the most significant limb. Use temporary storage not to clobber
92
the original contents of the denominator. */
93
tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
94
mpn_lshift (tp, dp, dsize, normalization_steps);
99
MPN_ZERO (rp, rsize - nsize);
100
nlimb = mpn_lshift (rp + (rsize - nsize),
101
np, nsize, normalization_steps);
105
nlimb = mpn_lshift (rp, np + (nsize - rsize),
106
rsize, normalization_steps);
138
/* separate blocks, for malloc debugging */
139
remp = TMP_ALLOC_LIMBS (dsize);
140
tp = (zeros > 0 ? TMP_ALLOC_LIMBS (tsize) : NULL);
118
MPN_ZERO (rp, rsize - nsize);
119
MPN_COPY (rp + (rsize - nsize), np, nsize);
123
MPN_COPY (rp, np + (nsize - rsize), rsize);
144
/* one block with conditionalized size, for efficiency */
145
remp = TMP_ALLOC_LIMBS (dsize + (zeros > 0 ? tsize : 0));
127
qlimb = mpn_divmod (qp, rp, rsize, dp, dsize);
128
qsize = rsize - dsize;
149
/* zero extend n into temporary space, if necessary */
152
MPN_ZERO (tp, zeros);
153
MPN_COPY (tp+zeros, np, nsize);
138
mp_size_t scale = nsize - dsize - N_QLIMBS;
140
#if defined (__vax__)
141
/* Ignore excess quotient limbs. This is necessary on a vax
142
with its small double exponent, since we'd otherwise get
143
exponent overflow while forming RES. */
144
if (qsize > N_QLIMBS)
146
qp += qsize - N_QLIMBS;
147
scale += qsize - N_QLIMBS;
153
for (i = qsize - 2; i >= 0; i--)
154
res = res * MP_BASE_AS_DOUBLE + qp[i];
156
res = __gmp_scale2 (res, GMP_NUMB_BITS * scale);
159
return sign_quotient >= 0 ? res : -res;
158
ASSERT (qsize == nsize - dsize + 1);
159
mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize);
161
/* strip possible zero high limb */
162
qsize -= (qp[qsize-1] == 0);
164
res = mpn_get_d (qp, qsize, sign_quotient, exp);
added
removed
src/gmp/mpq/get_d.c
