get_d.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:13k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpn_get_d -- limbs to double conversion.
  2.    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
  3.    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
  4.    FUTURE GNU MP RELEASES.
  5. Copyright 2003, 2004, 2007, 2009 Free Software Foundation, Inc.
  6. 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 Lesser General Public License as published by
  9. the Free Software Foundation; either version 3 of the License, or (at your
  10. option) any later version.
  11. The GNU MP Library is distributed in the hope that it will be useful, but
  12. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  13. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  14. License for more details.
  15. You should have received a copy of the GNU Lesser General Public License
  16. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  17. #include "gmp.h"
  18. #include "gmp-impl.h"
  19. #include "longlong.h"
  20. #ifndef _GMP_IEEE_FLOATS
  21. #define _GMP_IEEE_FLOATS 0
  22. #endif
  23. #if ! _GMP_IEEE_FLOATS
  24. /* dummy definition, just to let dead code compile */
  25. union ieee_double_extract {
  26.   struct {
  27.     int manh, manl, sig, exp;
  28.   } s;
  29.   double d;
  30. };
  31. #endif
  32. /* To force use of the generic C code for testing, put
  33.    "#define _GMP_IEEE_FLOATS 0" at this point.  */
  34. /* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
  35.    rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
  36.    wrong if that addition overflows.
  37.    The workaround here avoids this bug by ensuring n is not a literal
  38.    constant.  Note that this is alpha specific.  The offending transformation
  39.    is/was in alpha.c alpha_emit_conditional_branch() under "We want to use
  40.    cmpcc/bcc".
  41.    Bizarrely, it turns out this happens also with Cray cc on
  42.    alphaev5-cray-unicosmk2.0.6.X, and has the same solution.  Don't know why
  43.    or how.  */
  44. #if HAVE_HOST_CPU_FAMILY_alpha
  45.   && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4))
  46.       || defined (_CRAY))
  47. static volatile const long CONST_1024 = 1024;
  48. static volatile const long CONST_NEG_1023 = -1023;
  49. static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
  50. #else
  51. #define CONST_1024       (1024)
  52. #define CONST_NEG_1023       (-1023)
  53. #define CONST_NEG_1022_SUB_53 (-1022 - 53)
  54. #endif
  55. /* Return the value {ptr,size}*2^exp, and negative if sign<0.
  56.    Must have size>=1, and a non-zero high limb ptr[size-1].
  57.    {ptr,size} is truncated towards zero.  This is consistent with other gmp
  58.    conversions, like mpz_set_f or mpz_set_q, and is easy to implement and
  59.    test.
  60.    In the past conversions had attempted (imperfectly) to let the hardware
  61.    float rounding mode take effect, but that gets tricky since multiple
  62.    roundings need to be avoided, or taken into account, and denorms mean the
  63.    effective precision of the mantissa is not constant.  (For reference,
  64.    mpz_get_d on IEEE systems was ok, except it operated on the absolute
  65.    value.  mpf_get_d and mpq_get_d suffered from multiple roundings and from
  66.    not always using enough bits to get the rounding right.)
  67.    It's felt that GMP is not primarily concerned with hardware floats, and
  68.    really isn't enhanced by getting involved with hardware rounding modes
  69.    (which could even be some weird unknown style), so something unambiguous
  70.    and straightforward is best.
  71.    The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
  72.    limb and is done with shifts and masks.  The 64-bit case in particular
  73.    should come out nice and compact.
  74.    The generic code works one bit at a time, which will be quite slow, but
  75.    should support any binary-based "double" and be safe against any rounding
  76.    mode.  Note in particular it works on IEEE systems too.
  77.    Traps:
  78.    Hardware traps for overflow to infinity, underflow to zero, or
  79.    unsupported denorms may or may not be taken.  The IEEE code works bitwise
  80.    and so probably won't trigger them, the generic code works by float
  81.    operations and so probably will.  This difference might be thought less
  82.    than ideal, but again its felt straightforward code is better than trying
  83.    to get intimate with hardware exceptions (of perhaps unknown nature).
  84.    Not done:
  85.    mpz_get_d in the past handled size==1 with a cast limb->double.  This
  86.    might still be worthwhile there (for up to the mantissa many bits), but
  87.    for mpn_get_d here, the cost of applying "exp" to the resulting exponent
  88.    would probably use up any benefit a cast may have over bit twiddling.
  89.    Also, if the exponent is pushed into denorm range then bit twiddling is
  90.    the only option, to ensure the desired truncation is obtained.
  91.    Other:
  92.    For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
  93.    to the kernel for values >= 2^63.  This makes it slow, and worse the kernel
  94.    Linux (what versions?) apparently uses untested code in its trap handling
  95.    routines, and gets the sign wrong.  We don't use such a limb-to-double
  96.    cast, neither in the IEEE or generic code.  */
  97. double
  98. mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
  99. {
  100.   ASSERT (size >= 0);
  101.   ASSERT_MPN (up, size);
  102.   ASSERT (size == 0 || up[size-1] != 0);
  103.   if (size == 0)
  104.     return 0.0;
  105.   /* Adjust exp to a radix point just above {up,size}, guarding against
  106.      overflow.  After this exp can of course be reduced to anywhere within
  107.      the {up,size} region without underflow.  */
  108.   if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
  109. > (unsigned long) (LONG_MAX - exp)))
  110.     {
  111.       if (_GMP_IEEE_FLOATS)
  112. goto ieee_infinity;
  113.       /* generic */
  114.       exp = LONG_MAX;
  115.     }
  116.   else
  117.     {
  118.       exp += GMP_NUMB_BITS * size;
  119.     }
  120. #if 1
  121. {
  122.   int lshift, nbits;
  123.   union ieee_double_extract u;
  124.   mp_limb_t x, mhi, mlo;
  125. #if GMP_LIMB_BITS == 64
  126.   mp_limb_t m;
  127.   up += size;
  128.   m = *--up;
  129.   count_leading_zeros (lshift, m);
  130.   exp -= (lshift - GMP_NAIL_BITS) + 1;
  131.   m <<= lshift;
  132.   nbits = GMP_LIMB_BITS - lshift;
  133.   if (nbits < 53 && size > 1)
  134.     {
  135.       x = *--up;
  136.       x <<= GMP_NAIL_BITS;
  137.       x >>= nbits;
  138.       m |= x;
  139.       nbits += GMP_NUMB_BITS;
  140.       if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
  141. {
  142.   x = *--up;
  143.   x <<= GMP_NAIL_BITS;
  144.   x >>= nbits;
  145.   m |= x;
  146.   nbits += GMP_NUMB_BITS;
  147. }
  148.     }
  149.   mhi = m >> (32 + 11);
  150.   mlo = m >> 11;
  151. #endif
  152. #if GMP_LIMB_BITS == 32
  153.   up += size;
  154.   x = *--up, size--;
  155.   count_leading_zeros (lshift, x);
  156.   exp -= (lshift - GMP_NAIL_BITS) + 1;
  157.   x <<= lshift;
  158.   mhi = x >> 11;
  159.   if (lshift < 11) /* FIXME: never true if NUMB < 20 bits */
  160.     {
  161.       /* All 20 bits in mhi */
  162.       mlo = x << 21;
  163.       /* >= 1 bit in mlo */
  164.       nbits = GMP_LIMB_BITS - lshift - 21;
  165.     }
  166.   else
  167.     {
  168.       if (size != 0)
  169. {
  170.   nbits = GMP_LIMB_BITS - lshift;
  171.   x = *--up, size--;
  172.   x <<= GMP_NAIL_BITS;
  173.   mhi |= x >> nbits >> 11;
  174.   mlo = x << GMP_LIMB_BITS - nbits - 11;
  175.   nbits = nbits + 11 - GMP_NAIL_BITS;
  176. }
  177.       else
  178. {
  179.   mlo = 0;
  180.   goto done;
  181. }
  182.     }
  183.   if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size != 0)
  184.     {
  185.       x = *--up, size--;
  186.       x <<= GMP_NAIL_BITS;
  187.       x >>= nbits;
  188.       mlo |= x;
  189.       nbits += GMP_NUMB_BITS;
  190.       if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size != 0)
  191. {
  192.   x = *--up, size--;
  193.   x <<= GMP_NAIL_BITS;
  194.   x >>= nbits;
  195.   mlo |= x;
  196.   nbits += GMP_NUMB_BITS;
  197.   if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size != 0)
  198.     {
  199.       x = *--up;
  200.       x <<= GMP_NAIL_BITS;
  201.       x >>= nbits;
  202.       mlo |= x;
  203.       nbits += GMP_NUMB_BITS;
  204.     }
  205. }
  206.     }
  207.  done:;
  208. #endif
  209.   {
  210.     if (UNLIKELY (exp >= CONST_1024))
  211.       {
  212. /* overflow, return infinity */
  213.       ieee_infinity:
  214. mhi = 0;
  215. mlo = 0;
  216. exp = 1024;
  217.       }
  218.     else if (UNLIKELY (exp <= CONST_NEG_1023))
  219.       {
  220. int rshift;
  221. if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
  222.   return 0.0;  /* denorm underflows to zero */
  223. rshift = -1022 - exp;
  224. ASSERT (rshift > 0 && rshift < 53);
  225. #if GMP_LIMB_BITS > 53
  226. mlo >>= rshift;
  227. mhi = mlo >> 32;
  228. #else
  229. if (rshift >= 32)
  230.   {
  231.     mlo = mhi;
  232.     mhi = 0;
  233.     rshift -= 32;
  234.   }
  235. lshift = GMP_LIMB_BITS - rshift;
  236. mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
  237. mhi >>= rshift;
  238. #endif
  239. exp = -1023;
  240.       }
  241.   }
  242.   u.s.manh = mhi;
  243.   u.s.manl = mlo;
  244.   u.s.exp = exp + 1023;
  245.   u.s.sig = (sign < 0);
  246.   return u.d;
  247. }
  248. #else
  249. #define ONE_LIMB    (GMP_LIMB_BITS == 64 && 2*GMP_NUMB_BITS >= 53)
  250. #define TWO_LIMBS   (GMP_LIMB_BITS == 32 && 3*GMP_NUMB_BITS >= 53)
  251.   if (_GMP_IEEE_FLOATS && (ONE_LIMB || TWO_LIMBS))
  252.     {
  253.       union ieee_double_extract  u;
  254.       mp_limb_t  m0, m1, m2, rmask;
  255.       int  lshift, rshift;
  256.       m0 = up[size-1];     /* high limb */
  257.       m1 = (size >= 2 ? up[size-2] : 0);   /* second highest limb */
  258.       count_leading_zeros (lshift, m0);
  259.       /* relative to just under high non-zero bit */
  260.       exp -= (lshift - GMP_NAIL_BITS) + 1;
  261.       if (ONE_LIMB)
  262. {
  263.   /* lshift to have high of m0 non-zero, and collapse nails */
  264.   rshift = GMP_LIMB_BITS - lshift;
  265.   m1 <<= GMP_NAIL_BITS;
  266.   rmask = GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX;
  267.   m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
  268.   /* rshift back to have bit 53 of m0 the high non-zero */
  269.   m0 >>= 11;
  270. }
  271.       else /* TWO_LIMBS */
  272. {
  273.   m2 = (size >= 3 ? up[size-3] : 0);  /* third highest limb */
  274.   /* collapse nails from m1 and m2 */
  275. #if GMP_NAIL_BITS != 0
  276.   m1 = (m1 << GMP_NAIL_BITS) | (m2 >> (GMP_NUMB_BITS-GMP_NAIL_BITS));
  277.   m2 <<= 2*GMP_NAIL_BITS;
  278. #endif
  279.   /* lshift to have high of m0:m1 non-zero, collapse nails from m0 */
  280.   rshift = GMP_LIMB_BITS - lshift;
  281.   rmask = (GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX);
  282.   m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
  283.   m1 = (m1 << lshift) | ((m2 >> rshift) & rmask);
  284.   /* rshift back to have bit 53 of m0:m1 the high non-zero */
  285.   m1 = (m1 >> 11) | (m0 << (GMP_LIMB_BITS-11));
  286.   m0 >>= 11;
  287. }
  288.       if (UNLIKELY (exp >= CONST_1024))
  289. {
  290.   /* overflow, return infinity */
  291. ieee_infinity:
  292.   m0 = 0;
  293.   m1 = 0;
  294.   exp = 1024;
  295. }
  296.       else if (UNLIKELY (exp <= CONST_NEG_1023))
  297. {
  298.   if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
  299.     return 0.0;  /* denorm underflows to zero */
  300.   rshift = -1022 - exp;
  301.   ASSERT (rshift > 0 && rshift < 53);
  302.   if (ONE_LIMB)
  303.     {
  304.       m0 >>= rshift;
  305.     }
  306.   else /* TWO_LIMBS */
  307.     {
  308.       if (rshift >= 32)
  309. {
  310.   m1 = m0;
  311.   m0 = 0;
  312.   rshift -= 32;
  313. }
  314.       lshift = GMP_LIMB_BITS - rshift;
  315.       m1 = (m1 >> rshift) | (rshift == 0 ? 0 : m0 << lshift);
  316.       m0 >>= rshift;
  317.     }
  318.   exp = -1023;
  319. }
  320.       if (ONE_LIMB)
  321. {
  322. #if GMP_LIMB_BITS > 32 /* avoid compiler warning about big shift */
  323.   u.s.manh = m0 >> 32;
  324. #endif
  325.   u.s.manl = m0;
  326. }
  327.       else /* TWO_LIMBS */
  328. {
  329.   u.s.manh = m0;
  330.   u.s.manl = m1;
  331. }
  332.       u.s.exp = exp + 1023;
  333.       u.s.sig = (sign < 0);
  334.       return u.d;
  335.     }
  336.   else
  337.     {
  338.       /* Non-IEEE or strange limb size, do something generic. */
  339.       mp_size_t      i;
  340.       mp_limb_t      limb, bit;
  341.       int      shift;
  342.       double      base, factor, prev_factor, d, new_d, diff;
  343.       /* "limb" is "up[i]" the limb being examined, "bit" is a mask for the
  344.  bit being examined, initially the highest non-zero bit.  */
  345.       i = size-1;
  346.       limb = up[i];
  347.       count_leading_zeros (shift, limb);
  348.       bit = GMP_LIMB_HIGHBIT >> shift;
  349.       /* relative to just under high non-zero bit */
  350.       exp -= (shift - GMP_NAIL_BITS) + 1;
  351.       /* Power up "factor" to 2^exp, being the value of the "bit" in "limb"
  352.  being examined.  */
  353.       base = (exp >= 0 ? 2.0 : 0.5);
  354.       exp = ABS (exp);
  355.       factor = 1.0;
  356.       for (;;)
  357. {
  358.   if (exp & 1)
  359.     {
  360.       prev_factor = factor;
  361.       factor *= base;
  362.       FORCE_DOUBLE (factor);
  363.       if (factor == 0.0)
  364. return 0.0; /* underflow */
  365.       if (factor == prev_factor)
  366. {
  367.   d = factor;   /* overflow, apparent infinity */
  368.   goto generic_done;
  369. }
  370.     }
  371.   exp >>= 1;
  372.   if (exp == 0)
  373.     break;
  374.   base *= base;
  375. }
  376.       /* Add a "factor" for each non-zero bit, working from high to low.
  377.  Stop if any rounding occurs, hence implementing a truncation.
  378.  Note no attention is paid to DBL_MANT_DIG, since the effective
  379.  number of bits in the mantissa isn't constant when in denorm range.
  380.  We also encountered an ARM system with apparently somewhat doubtful
  381.  software floats where DBL_MANT_DIG claimed 53 bits but only 32
  382.  actually worked.  */
  383.       d = factor;  /* high bit */
  384.       for (;;)
  385. {
  386.   factor *= 0.5;  /* next bit */
  387.   bit >>= 1;
  388.   if (bit == 0)
  389.     {
  390.       /* next limb, if any */
  391.       i--;
  392.       if (i < 0)
  393. break;
  394.       limb = up[i];
  395.       bit = GMP_NUMB_HIGHBIT;
  396.     }
  397.   if (bit & limb)
  398.     {
  399.       new_d = d + factor;
  400.       FORCE_DOUBLE (new_d);
  401.       diff = new_d - d;
  402.       if (diff != factor)
  403. break;  /* rounding occured, stop now */
  404.       d = new_d;
  405.     }
  406. }
  407.     generic_done:
  408.       return (sign >= 0 ? d : -d);
  409.     }
  410. #endif
  411. }