1 /* mpf_set_q (mpf_t rop, mpq_t op) -- Convert the rational op to the float rop.
3 Copyright (C) 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 mpf_set_q (mpf_t r, mpq_srcptr q)
62 mp_size_t nsize, dsize;
63 mp_size_t qsize, rsize;
64 mp_size_t sign_quotient;
65 unsigned normalization_steps;
72 nsize = SIZ (&q->_mp_num);
73 dsize = SIZ (&q->_mp_den);
88 sign_quotient = nsize ^ dsize;
91 np = PTR (&q->_mp_num);
92 dp = PTR (&q->_mp_den);
107 rsize = MAX (nsize, dsize);
108 rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB);
110 count_leading_zeros (normalization_steps, dp[dsize - 1]);
112 /* Normalize the denominator, i.e. make its most significant bit set by
113 shifting it NORMALIZATION_STEPS bits to the left. Also shift the
114 numerator the same number of steps (to keep the quotient the same!). */
115 if (normalization_steps != 0)
120 /* Shift up the denominator setting the most significant bit of
121 the most significant limb. Use temporary storage not to clobber
122 the original contents of the denominator. */
123 tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
124 mpn_lshift (tp, dp, dsize, normalization_steps);
129 MPN_ZERO (rp, rsize - nsize);
130 nlimb = mpn_lshift (rp + (rsize - nsize),
131 np, nsize, normalization_steps);
135 nlimb = mpn_lshift (rp, np, nsize, normalization_steps);
148 MPN_ZERO (rp, rsize - nsize);
149 MPN_COPY (rp + (rsize - nsize), np, nsize);
153 MPN_COPY (rp, np, rsize);
157 qlimb = mpn_divrem (qp, prec - 1 - (rsize - dsize), rp, rsize, dp, dsize);