1 /* double mpq_get_d (mpq_t src) -- Return the double approximation to SRC.
3 Copyright (C) 1995, 1996 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 Library General Public License as published by
9 the Free Software Foundation; either version 2 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 Library General Public
15 License for more details.
17 You should have received a copy of the GNU Library 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 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:
35 0. Let n be somewhat larger than the number of significant bits in a double.
36 1. Extract the most significant n bits of the denominator, and an equal
37 number of bits from the numerator.
38 2. Interpret the extracted numbers as integers, call them a and b
39 respectively, and develop n bits of the fractions ((a + 1) / b) and
40 (a / (b + 1)) using mpn_divrem.
41 3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT,
42 we are done. If they are different, repeat the algorithm from step 1,
43 but first let n = n * 2.
44 4. If we end up using all bits from the numerator and denominator, fall
45 back to the first algorithm above.
46 5. Just to make life harder, The computation of a + 1 and b + 1 above
47 might give carry-out... Needs special handling. It might work to
48 subtract 1 in both cases instead.
53 mpq_get_d (const MP_RAT *src)
61 mp_size_t nsize = src->_mp_num._mp_size;
62 mp_size_t dsize = src->_mp_den._mp_size;
63 mp_size_t qsize, rsize;
64 mp_size_t sign_quotient = nsize ^ dsize;
65 unsigned normalization_steps;
67 #define N_QLIMBS (1 + (sizeof (double) + BYTES_PER_MP_LIMB-1) / BYTES_PER_MP_LIMB)
68 mp_limb_t qp[N_QLIMBS + 1];
77 np = src->_mp_num._mp_d;
78 dp = src->_mp_den._mp_d;
80 rsize = dsize + N_QLIMBS;
81 rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB);
83 count_leading_zeros (normalization_steps, dp[dsize - 1]);
85 /* Normalize the denominator, i.e. make its most significant bit set by
86 shifting it NORMALIZATION_STEPS bits to the left. Also shift the
87 numerator the same number of steps (to keep the quotient the same!). */
88 if (normalization_steps != 0)
93 /* Shift up the denominator setting the most significant bit of
94 the most significant limb. Use temporary storage not to clobber
95 the original contents of the denominator. */
96 tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
97 mpn_lshift (tp, dp, dsize, normalization_steps);
102 MPN_ZERO (rp, rsize - nsize);
103 nlimb = mpn_lshift (rp + (rsize - nsize),
104 np, nsize, normalization_steps);
108 nlimb = mpn_lshift (rp, np + (nsize - rsize),
109 rsize, normalization_steps);
121 MPN_ZERO (rp, rsize - nsize);
122 MPN_COPY (rp + (rsize - nsize), np, nsize);
126 MPN_COPY (rp, np + (nsize - rsize), rsize);
130 qlimb = mpn_divmod (qp, rp, rsize, dp, dsize);
131 qsize = rsize - dsize;
143 for (i = qsize - 2; i >= 0; i--)
144 res = res * MP_BASE_AS_DOUBLE + qp[i];
146 res = __gmp_scale2 (res, BITS_PER_MP_LIMB * (nsize - dsize - N_QLIMBS));
149 return sign_quotient >= 0 ? res : -res;