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