diff options
Diffstat (limited to 'vendor/gmp-6.3.0/mpq/get_d.c')
| -rw-r--r-- | vendor/gmp-6.3.0/mpq/get_d.c | 162 |
1 files changed, 162 insertions, 0 deletions
diff --git a/vendor/gmp-6.3.0/mpq/get_d.c b/vendor/gmp-6.3.0/mpq/get_d.c new file mode 100644 index 0000000..4d9779d --- /dev/null +++ b/vendor/gmp-6.3.0/mpq/get_d.c @@ -0,0 +1,162 @@ +/* double mpq_get_d (mpq_t src) -- mpq to double, rounding towards zero. + +Copyright 1995, 1996, 2001-2005, 2018, 2019 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> /* for NULL */ +#include "gmp-impl.h" +#include "longlong.h" + + +/* All that's needed is to get the high 53 bits of the quotient num/den, + rounded towards zero. More than 53 bits is fine, any excess is ignored + by mpn_get_d. + + N_QLIMBS is how many quotient limbs we need to satisfy the mantissa of a + double, assuming the highest of those limbs is non-zero. The target + qsize for mpn_tdiv_qr is then 1 more than this, since that function may + give a zero in the high limb (and non-zero in the second highest). + + The use of 8*sizeof(double) in N_QLIMBS is an overestimate of the + mantissa bits, but it gets the same result as the true value (53 or 48 or + whatever) when rounded up to a multiple of GMP_NUMB_BITS, for non-nails. + + Enhancements: + + Use the true mantissa size in the N_QLIMBS formula, to save a divide step + in nails. + + Examine the high limbs of num and den to see if the highest 1 bit of the + quotient will fall high enough that just N_QLIMBS-1 limbs is enough to + get the necessary bits, thereby saving a division step. + + Bit shift either num or den to arrange for the above condition on the + high 1 bit of the quotient, to save a division step always. A shift to + save a division step is definitely worthwhile with mpn_tdiv_qr, though we + may want to reassess this on big num/den when a quotient-only division + exists. + + Maybe we could estimate the final exponent using nsize-dsize (and + possibly the high limbs of num and den), so as to detect overflow and + return infinity or zero quickly. Overflow is never very helpful to an + application, and can therefore probably be regarded as abnormal, but we + may still like to optimize it if the conditions are easy. (This would + only be for float formats we know, unknown formats are not important and + can be left to mpn_get_d.) + + Future: + + If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of + padding n with zeros in temporary space. + + Alternatives: + + An alternative algorithm, that may be faster: + 0. Let n be somewhat larger than the number of significant bits in a double. + 1. Extract the most significant n bits of the denominator, and an equal + number of bits from the numerator. + 2. Interpret the extracted numbers as integers, call them a and b + respectively, and develop n bits of the fractions ((a + 1) / b) and + (a / (b + 1)) using mpn_divrem. + 3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT, + we are done. If they are different, repeat the algorithm from step 1, + but first let n = n * 2. + 4. If we end up using all bits from the numerator and denominator, fall + back to a plain division. + 5. Just to make life harder, The computation of a + 1 and b + 1 above + might give carry-out... Needs special handling. It might work to + subtract 1 in both cases instead. + + Not certain if this approach would be faster than a quotient-only + division. Presumably such optimizations are the sort of thing we would + like to have helping everywhere that uses a quotient-only division. */ + +double +mpq_get_d (mpq_srcptr src) +{ + double res; + mp_srcptr np, dp; + mp_ptr temp; + mp_size_t nsize = SIZ(NUM(src)); + mp_size_t dsize = SIZ(DEN(src)); + mp_size_t qsize, prospective_qsize, zeros; + mp_size_t sign_quotient = nsize; + long exp; +#define N_QLIMBS (1 + (sizeof (double) + GMP_LIMB_BYTES-1) / GMP_LIMB_BYTES) + mp_limb_t qarr[N_QLIMBS + 1]; + mp_ptr qp = qarr; + TMP_DECL; + + ASSERT (dsize > 0); /* canonical src */ + + /* mpn_get_d below requires a non-zero operand */ + if (UNLIKELY (nsize == 0)) + return 0.0; + + TMP_MARK; + nsize = ABS (nsize); + dsize = ABS (dsize); + np = PTR(NUM(src)); + dp = PTR(DEN(src)); + + prospective_qsize = nsize - dsize; /* from using given n,d */ + qsize = N_QLIMBS; /* desired qsize */ + + zeros = qsize - prospective_qsize; /* padding n to get qsize */ + exp = (long) -zeros * GMP_NUMB_BITS; /* relative to low of qp */ + + /* zero extend n into temporary space, if necessary */ + if (zeros > 0) + { + mp_size_t tsize; + tsize = nsize + zeros; /* size for copy of n */ + + temp = TMP_ALLOC_LIMBS (tsize + 1); + MPN_FILL (temp, zeros, 0); + MPN_COPY (temp + zeros, np, nsize); + np = temp; + nsize = tsize; + } + else /* negative zeros means shorten n */ + { + np -= zeros; + nsize += zeros; + + temp = TMP_ALLOC_LIMBS (nsize + 1); + } + + ASSERT (qsize == nsize - dsize); + mpn_div_q (qp, np, nsize, dp, dsize, temp); + + /* strip possible zero high limb */ + qsize += (qp[qsize] != 0); + + res = mpn_get_d (qp, qsize, sign_quotient, exp); + TMP_FREE; + return res; +} |