diff options
Diffstat (limited to 'vendor/gmp-6.3.0/mpn/generic/get_d.c')
-rw-r--r-- | vendor/gmp-6.3.0/mpn/generic/get_d.c | 438 |
1 files changed, 438 insertions, 0 deletions
diff --git a/vendor/gmp-6.3.0/mpn/generic/get_d.c b/vendor/gmp-6.3.0/mpn/generic/get_d.c new file mode 100644 index 0000000..8bef128 --- /dev/null +++ b/vendor/gmp-6.3.0/mpn/generic/get_d.c @@ -0,0 +1,438 @@ +/* mpn_get_d -- limbs to double conversion. + + THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST + CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN + FUTURE GNU MP RELEASES. + +Copyright 2003, 2004, 2007, 2009, 2010, 2012, 2018 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 "config.h" + +#if HAVE_FLOAT_H +#include <float.h> /* for DBL_MANT_DIG and FLT_RADIX */ +#endif + +#include "gmp-impl.h" +#include "longlong.h" + +#ifndef _GMP_IEEE_FLOATS +#define _GMP_IEEE_FLOATS 0 +#endif + +/* To force use of the generic C code for testing, put + "#define _GMP_IEEE_FLOATS 0" at this point. */ + + +/* In alpha gcc prior to 3.4, signed DI comparisons involving constants are + rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly + wrong if that addition overflows. + + The workaround here avoids this bug by ensuring n is not a literal constant. + Note that this is alpha specific. The offending transformation is/was in + alpha.c alpha_emit_conditional_branch() under "We want to use cmpcc/bcc". + + Bizarrely, this happens also with Cray cc on alphaev5-cray-unicosmk2.0.6.X, + and has the same solution. Don't know why or how. */ + +#if HAVE_HOST_CPU_FAMILY_alpha \ + && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4)) \ + || defined (_CRAY)) +static volatile const long CONST_1024 = 1024; +static volatile const long CONST_NEG_1023 = -1023; +static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53; +#else +#define CONST_1024 (1024) +#define CONST_NEG_1023 (-1023) +#define CONST_NEG_1022_SUB_53 (-1022 - 53) +#endif + + +/* Return the value {ptr,size}*2^exp, and negative if sign<0. Must have + size>=1, and a non-zero high limb ptr[size-1]. + + When we know the fp format, the result is truncated towards zero. This is + consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is + easy to implement and test. + + When we do not know the format, such truncation seems much harder. One + would need to defeat any rounding mode, including round-up. + + It's felt that GMP is not primarily concerned with hardware floats, and + really isn't enhanced by getting involved with hardware rounding modes + (which could even be some weird unknown style), so something unambiguous and + straightforward is best. + + + The IEEE code below is the usual case, it knows either a 32-bit or 64-bit + limb and is done with shifts and masks. The 64-bit case in particular + should come out nice and compact. + + The generic code used to work one bit at a time, which was not only slow, + but implicitly relied upon denorms for intermediates, since the lowest bits' + weight of a perfectly valid fp number underflows in non-denorm. Therefore, + the generic code now works limb-per-limb, initially creating a number x such + that 1 <= x <= BASE. (BASE is reached only as result of rounding.) Then + x's exponent is scaled with explicit code (not ldexp to avoid libm + dependency). It is a tap-dance to avoid underflow or overflow, beware! + + + Traps: + + Hardware traps for overflow to infinity, underflow to zero, or unsupported + denorms may or may not be taken. The IEEE code works bitwise and so + probably won't trigger them, the generic code works by float operations and + so probably will. This difference might be thought less than ideal, but + again its felt straightforward code is better than trying to get intimate + with hardware exceptions (of perhaps unknown nature). + + + Not done: + + mpz_get_d in the past handled size==1 with a cast limb->double. This might + still be worthwhile there (for up to the mantissa many bits), but for + mpn_get_d here, the cost of applying "exp" to the resulting exponent would + probably use up any benefit a cast may have over bit twiddling. Also, if + the exponent is pushed into denorm range then bit twiddling is the only + option, to ensure the desired truncation is obtained. + + + Other: + + For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL + to the kernel for values >= 2^63. This makes it slow, and worse the kernel + Linux (what versions?) apparently uses untested code in its trap handling + routines, and gets the sign wrong. We don't use such a limb-to-double + cast, neither in the IEEE or generic code. */ + + + +#undef FORMAT_RECOGNIZED + +double +mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp) +{ + int lshift, nbits; + mp_limb_t x, mhi, mlo; + + ASSERT (size >= 0); + ASSERT_MPN (up, size); + ASSERT (size == 0 || up[size-1] != 0); + + if (size == 0) + return 0.0; + + /* Adjust exp to a radix point just above {up,size}, guarding against + overflow. After this exp can of course be reduced to anywhere within + the {up,size} region without underflow. */ + if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size) + > ((unsigned long) LONG_MAX - exp))) + { +#if _GMP_IEEE_FLOATS + goto ieee_infinity; +#endif + + /* generic */ + exp = LONG_MAX; + } + else + { + exp += GMP_NUMB_BITS * size; + } + +#if _GMP_IEEE_FLOATS + { + union ieee_double_extract u; + + up += size; + +#if GMP_LIMB_BITS == 64 + mlo = up[-1]; + count_leading_zeros (lshift, mlo); + + exp -= (lshift - GMP_NAIL_BITS) + 1; + mlo <<= lshift; + + nbits = GMP_LIMB_BITS - lshift; + + if (nbits < 53 && size > 1) + { + x = up[-2]; + x <<= GMP_NAIL_BITS; + x >>= nbits; + mlo |= x; + nbits += GMP_NUMB_BITS; + + if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2) + { + x = up[-3]; + x <<= GMP_NAIL_BITS; + x >>= nbits; + mlo |= x; + nbits += GMP_NUMB_BITS; + } + } + mhi = mlo >> (32 + 11); + mlo = mlo >> 11; /* later implicitly truncated to 32 bits */ +#endif +#if GMP_LIMB_BITS == 32 + x = *--up; + count_leading_zeros (lshift, x); + + exp -= (lshift - GMP_NAIL_BITS) + 1; + x <<= lshift; + mhi = x >> 11; + + if (lshift < 11) /* FIXME: never true if NUMB < 20 bits */ + { + /* All 20 bits in mhi */ + mlo = x << 21; + /* >= 1 bit in mlo */ + nbits = GMP_LIMB_BITS - lshift - 21; + } + else + { + if (size > 1) + { + nbits = GMP_LIMB_BITS - lshift; + + x = *--up, size--; + x <<= GMP_NAIL_BITS; + mhi |= x >> nbits >> 11; + + mlo = x << (GMP_LIMB_BITS - nbits - 11); + nbits = nbits + 11 - GMP_NAIL_BITS; + } + else + { + mlo = 0; + goto done; + } + } + + /* Now all needed bits in mhi have been accumulated. Add bits to mlo. */ + + if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size > 1) + { + x = up[-1]; + x <<= GMP_NAIL_BITS; + x >>= nbits; + mlo |= x; + nbits += GMP_NUMB_BITS; + + if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size > 2) + { + x = up[-2]; + x <<= GMP_NAIL_BITS; + x >>= nbits; + mlo |= x; + nbits += GMP_NUMB_BITS; + + if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size > 3) + { + x = up[-3]; + x <<= GMP_NAIL_BITS; + x >>= nbits; + mlo |= x; + nbits += GMP_NUMB_BITS; + } + } + } + + done:; + +#endif + if (UNLIKELY (exp >= CONST_1024)) + { + /* overflow, return infinity */ + ieee_infinity: + mhi = 0; + mlo = 0; + exp = 1024; + } + else if (UNLIKELY (exp <= CONST_NEG_1023)) + { + int rshift; + + if (LIKELY (exp <= CONST_NEG_1022_SUB_53)) + return 0.0; /* denorm underflows to zero */ + + rshift = -1022 - exp; + ASSERT (rshift > 0 && rshift < 53); +#if GMP_LIMB_BITS > 53 + mlo >>= rshift; + mhi = mlo >> 32; +#else + if (rshift >= 32) + { + mlo = mhi; + mhi = 0; + rshift -= 32; + } + lshift = GMP_LIMB_BITS - rshift; + mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift); + mhi >>= rshift; +#endif + exp = -1023; + } + u.s.manh = mhi; + u.s.manl = mlo; + u.s.exp = exp + 1023; + u.s.sig = (sign < 0); + return u.d; + } +#define FORMAT_RECOGNIZED 1 +#endif + +#if HAVE_DOUBLE_VAX_D + { + union double_extract u; + + up += size; + + mhi = up[-1]; + + count_leading_zeros (lshift, mhi); + exp -= lshift; + mhi <<= lshift; + + mlo = 0; + if (size > 1) + { + mlo = up[-2]; + if (lshift != 0) + mhi += mlo >> (GMP_LIMB_BITS - lshift); + mlo <<= lshift; + + if (size > 2 && lshift > 8) + { + x = up[-3]; + mlo += x >> (GMP_LIMB_BITS - lshift); + } + } + + if (UNLIKELY (exp >= 128)) + { + /* overflow, return maximum number */ + mhi = 0xffffffff; + mlo = 0xffffffff; + exp = 127; + } + else if (UNLIKELY (exp < -128)) + { + return 0.0; /* underflows to zero */ + } + + u.s.man3 = mhi >> 24; /* drop msb, since implicit */ + u.s.man2 = mhi >> 8; + u.s.man1 = (mhi << 8) + (mlo >> 24); + u.s.man0 = mlo >> 8; + u.s.exp = exp + 128; + u.s.sig = sign < 0; + return u.d; + } +#define FORMAT_RECOGNIZED 1 +#endif + +#if ! FORMAT_RECOGNIZED + +#if !defined(GMP_DBL_MANT_BITS) +#if defined(DBL_MANT_DIG) && FLT_RADIX == 2 +#define GMP_DBL_MANT_BITS DBL_MANT_DIG +#else +/* FIXME: Chose a smarter default value. */ +#define GMP_DBL_MANT_BITS (16 * sizeof (double)) +#endif +#endif + + { /* Non-IEEE or strange limb size, generically convert + GMP_DBL_MANT_BITS bits. */ + mp_limb_t l; + int m; + mp_size_t i; + double d, weight; + unsigned long uexp; + + /* First generate an fp number disregarding exp, instead keeping things + within the numb base factor from 1, which should prevent overflow and + underflow even for the most exponent limited fp formats. */ + i = size - 1; + l = up[i]; + count_leading_zeros (m, l); + m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS; + if (m < 0) + l &= GMP_NUMB_MAX << -m; + d = l; + for (weight = 1/MP_BASE_AS_DOUBLE; m > 0 && --i >= 0;) + { + l = up[i]; + m -= GMP_NUMB_BITS; + if (m < 0) + l &= GMP_NUMB_MAX << -m; + d += l * weight; + weight /= MP_BASE_AS_DOUBLE; + if (weight == 0) + break; + } + + /* Now apply exp. */ + exp -= GMP_NUMB_BITS; + if (exp > 0) + { + weight = 2.0; + uexp = exp; + } + else + { + weight = 0.5; + uexp = NEG_CAST (unsigned long, exp); + } +#if 1 + /* Square-and-multiply exponentiation. */ + if (uexp & 1) + d *= weight; + while (uexp >>= 1) + { + weight *= weight; + if (uexp & 1) + d *= weight; + } +#else + /* Plain exponentiation. */ + while (uexp > 0) + { + d *= weight; + uexp--; + } +#endif + + return sign >= 0 ? d : -d; + } +#endif +} |