1 /* double mpq_get_d (mpq_t src) -- mpq to double, rounding towards zero.
2
3 Copyright 1995, 1996, 2001-2005, 2018, 2019 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library.
6
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of either:
9
10 * the GNU Lesser General Public License as published by the Free
11 Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 or
15
16 * the GNU General Public License as published by the Free Software
17 Foundation; either version 2 of the License, or (at your option) any
18 later version.
19
20 or both in parallel, as here.
21
22 The GNU MP Library is distributed in the hope that it will be useful, but
23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25 for more details.
26
27 You should have received copies of the GNU General Public License and the
28 GNU Lesser General Public License along with the GNU MP Library. If not,
29 see https://www.gnu.org/licenses/. */
30
31 #include <stdio.h> /* for NULL */
32 #include "gmp-impl.h"
33 #include "longlong.h"
34
35
36 /* All that's needed is to get the high 53 bits of the quotient num/den,
37 rounded towards zero. More than 53 bits is fine, any excess is ignored
38 by mpn_get_d.
39
40 N_QLIMBS is how many quotient limbs we need to satisfy the mantissa of a
41 double, assuming the highest of those limbs is non-zero. The target
42 qsize for mpn_tdiv_qr is then 1 more than this, since that function may
43 give a zero in the high limb (and non-zero in the second highest).
44
45 The use of 8*sizeof(double) in N_QLIMBS is an overestimate of the
46 mantissa bits, but it gets the same result as the true value (53 or 48 or
47 whatever) when rounded up to a multiple of GMP_NUMB_BITS, for non-nails.
48
49 Enhancements:
50
51 Use the true mantissa size in the N_QLIMBS formula, to save a divide step
52 in nails.
53
54 Examine the high limbs of num and den to see if the highest 1 bit of the
55 quotient will fall high enough that just N_QLIMBS-1 limbs is enough to
56 get the necessary bits, thereby saving a division step.
57
58 Bit shift either num or den to arrange for the above condition on the
59 high 1 bit of the quotient, to save a division step always. A shift to
60 save a division step is definitely worthwhile with mpn_tdiv_qr, though we
61 may want to reassess this on big num/den when a quotient-only division
62 exists.
63
64 Maybe we could estimate the final exponent using nsize-dsize (and
65 possibly the high limbs of num and den), so as to detect overflow and
66 return infinity or zero quickly. Overflow is never very helpful to an
67 application, and can therefore probably be regarded as abnormal, but we
68 may still like to optimize it if the conditions are easy. (This would
69 only be for float formats we know, unknown formats are not important and
70 can be left to mpn_get_d.)
71
72 Future:
73
74 If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of
75 padding n with zeros in temporary space.
76
77 Alternatives:
78
79 An alternative algorithm, that may be faster:
80 0. Let n be somewhat larger than the number of significant bits in a double.
81 1. Extract the most significant n bits of the denominator, and an equal
82 number of bits from the numerator.
83 2. Interpret the extracted numbers as integers, call them a and b
84 respectively, and develop n bits of the fractions ((a + 1) / b) and
85 (a / (b + 1)) using mpn_divrem.
86 3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT,
87 we are done. If they are different, repeat the algorithm from step 1,
88 but first let n = n * 2.
89 4. If we end up using all bits from the numerator and denominator, fall
90 back to a plain division.
91 5. Just to make life harder, The computation of a + 1 and b + 1 above
92 might give carry-out... Needs special handling. It might work to
93 subtract 1 in both cases instead.
94
95 Not certain if this approach would be faster than a quotient-only
96 division. Presumably such optimizations are the sort of thing we would
97 like to have helping everywhere that uses a quotient-only division. */
98
99 double
100 mpq_get_d (mpq_srcptr src)
101 {
102 double res;
103 mp_srcptr np, dp;
104 mp_ptr temp;
105 mp_size_t nsize = SIZ(NUM(src));
106 mp_size_t dsize = SIZ(DEN(src));
107 mp_size_t qsize, prospective_qsize, zeros;
108 mp_size_t sign_quotient = nsize;
109 long exp;
110 #define N_QLIMBS (1 + (sizeof (double) + GMP_LIMB_BYTES-1) / GMP_LIMB_BYTES)
111 mp_limb_t qarr[N_QLIMBS + 1];
112 mp_ptr qp = qarr;
113 TMP_DECL;
114
115 ASSERT (dsize > 0); /* canonical src */
116
117 /* mpn_get_d below requires a non-zero operand */
118 if (UNLIKELY (nsize == 0))
119 return 0.0;
120
121 TMP_MARK;
122 nsize = ABS (nsize);
123 dsize = ABS (dsize);
124 np = PTR(NUM(src));
125 dp = PTR(DEN(src));
126
127 prospective_qsize = nsize - dsize; /* from using given n,d */
128 qsize = N_QLIMBS; /* desired qsize */
129
130 zeros = qsize - prospective_qsize; /* padding n to get qsize */
131 exp = (long) -zeros * GMP_NUMB_BITS; /* relative to low of qp */
132
133 /* zero extend n into temporary space, if necessary */
134 if (zeros > 0)
135 {
136 mp_size_t tsize;
137 tsize = nsize + zeros; /* size for copy of n */
138
139 temp = TMP_ALLOC_LIMBS (tsize + 1);
140 MPN_FILL (temp, zeros, 0);
141 MPN_COPY (temp + zeros, np, nsize);
142 np = temp;
143 nsize = tsize;
144 }
145 else /* negative zeros means shorten n */
146 {
147 np -= zeros;
148 nsize += zeros;
149
150 temp = TMP_ALLOC_LIMBS (nsize + 1);
151 }
152
153 ASSERT (qsize == nsize - dsize);
154 mpn_div_q (qp, np, nsize, dp, dsize, temp);
155
156 /* strip possible zero high limb */
157 qsize += (qp[qsize] != 0);
158
159 res = mpn_get_d (qp, qsize, sign_quotient, exp);
160 TMP_FREE;
161 return res;
162 }