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

数学计算

开发平台:

Unix_Linux

  1. /* __gmp_extract_double -- convert from double to array of mp_limb_t.
  2. Copyright 1996, 1999, 2000, 2001, 2002, 2006 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include "gmp.h"
  15. #include "gmp-impl.h"
  16. #ifdef XDEBUG
  17. #undef _GMP_IEEE_FLOATS
  18. #endif
  19. #ifndef _GMP_IEEE_FLOATS
  20. #define _GMP_IEEE_FLOATS 0
  21. #endif
  22. #define BITS_IN_MANTISSA 53
  23. /* Extract a non-negative double in d.  */
  24. int
  25. __gmp_extract_double (mp_ptr rp, double d)
  26. {
  27.   long exp;
  28.   unsigned sc;
  29. #ifdef _LONG_LONG_LIMB
  30. #define BITS_PER_PART 64 /* somewhat bogus */
  31.   unsigned long long int manl;
  32. #else
  33. #define BITS_PER_PART GMP_LIMB_BITS
  34.   unsigned long int manh, manl;
  35. #endif
  36.   /* BUGS
  37.      1. Should handle Inf and NaN in IEEE specific code.
  38.      2. Handle Inf and NaN also in default code, to avoid hangs.
  39.      3. Generalize to handle all GMP_LIMB_BITS >= 32.
  40.      4. This lits is incomplete and misspelled.
  41.    */
  42.   ASSERT (d >= 0.0);
  43.   if (d == 0.0)
  44.     {
  45.       MPN_ZERO (rp, LIMBS_PER_DOUBLE);
  46.       return 0;
  47.     }
  48. #if _GMP_IEEE_FLOATS
  49.   {
  50. #if defined (__alpha) && __GNUC__ == 2 && __GNUC_MINOR__ == 8
  51.     /* Work around alpha-specific bug in GCC 2.8.x.  */
  52.     volatile
  53. #endif
  54.     union ieee_double_extract x;
  55.     x.d = d;
  56.     exp = x.s.exp;
  57. #if BITS_PER_PART == 64 /* generalize this to BITS_PER_PART > BITS_IN_MANTISSA */
  58.     manl = (((mp_limb_t) 1 << 63)
  59.     | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
  60.     if (exp == 0)
  61.       {
  62. /* Denormalized number.  Don't try to be clever about this,
  63.    since it is not an important case to make fast.  */
  64. exp = 1;
  65. do
  66.   {
  67.     manl = manl << 1;
  68.     exp--;
  69.   }
  70. while ((manl & GMP_LIMB_HIGHBIT) == 0);
  71.       }
  72. #endif
  73. #if BITS_PER_PART == 32
  74.     manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
  75.     manl = x.s.manl << 11;
  76.     if (exp == 0)
  77.       {
  78. /* Denormalized number.  Don't try to be clever about this,
  79.    since it is not an important case to make fast.  */
  80. exp = 1;
  81. do
  82.   {
  83.     manh = (manh << 1) | (manl >> 31);
  84.     manl = manl << 1;
  85.     exp--;
  86.   }
  87. while ((manh & GMP_LIMB_HIGHBIT) == 0);
  88.       }
  89. #endif
  90. #if BITS_PER_PART != 32 && BITS_PER_PART != 64
  91.   You need to generalize the code above to handle this.
  92. #endif
  93.     exp -= 1022; /* Remove IEEE bias.  */
  94.   }
  95. #else
  96.   {
  97.     /* Unknown (or known to be non-IEEE) double format.  */
  98.     exp = 0;
  99.     if (d >= 1.0)
  100.       {
  101. ASSERT_ALWAYS (d * 0.5 != d);
  102. while (d >= 32768.0)
  103.   {
  104.     d *= (1.0 / 65536.0);
  105.     exp += 16;
  106.   }
  107. while (d >= 1.0)
  108.   {
  109.     d *= 0.5;
  110.     exp += 1;
  111.   }
  112.       }
  113.     else if (d < 0.5)
  114.       {
  115. while (d < (1.0 / 65536.0))
  116.   {
  117.     d *=  65536.0;
  118.     exp -= 16;
  119.   }
  120. while (d < 0.5)
  121.   {
  122.     d *= 2.0;
  123.     exp -= 1;
  124.   }
  125.       }
  126.     d *= (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
  127. #if BITS_PER_PART == 64
  128.     manl = d;
  129. #endif
  130. #if BITS_PER_PART == 32
  131.     manh = d;
  132.     manl = (d - manh) * (4.0 * ((unsigned long int) 1 << (BITS_PER_PART - 2)));
  133. #endif
  134.   }
  135. #endif /* IEEE */
  136.   sc = (unsigned) (exp + 64 * GMP_NUMB_BITS) % GMP_NUMB_BITS;
  137.   /* We add something here to get rounding right.  */
  138.   exp = (exp + 64 * GMP_NUMB_BITS) / GMP_NUMB_BITS - 64 * GMP_NUMB_BITS / GMP_NUMB_BITS + 1;
  139. #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 2
  140. #if GMP_NAIL_BITS == 0
  141.   if (sc != 0)
  142.     {
  143.       rp[1] = manl >> (GMP_LIMB_BITS - sc);
  144.       rp[0] = manl << sc;
  145.     }
  146.   else
  147.     {
  148.       rp[1] = manl;
  149.       rp[0] = 0;
  150.       exp--;
  151.     }
  152. #else
  153.   if (sc > GMP_NAIL_BITS)
  154.     {
  155.       rp[1] = manl >> (GMP_LIMB_BITS - sc);
  156.       rp[0] = (manl << (sc - GMP_NAIL_BITS)) & GMP_NUMB_MASK;
  157.     }
  158.   else
  159.     {
  160.       if (sc == 0)
  161. {
  162.   rp[1] = manl >> GMP_NAIL_BITS;
  163.   rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
  164.   exp--;
  165. }
  166.       else
  167. {
  168.   rp[1] = manl >> (GMP_LIMB_BITS - sc);
  169.   rp[0] = (manl >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
  170. }
  171.     }
  172. #endif
  173. #endif
  174. #if BITS_PER_PART == 64 && LIMBS_PER_DOUBLE == 3
  175.   if (sc > GMP_NAIL_BITS)
  176.     {
  177.       rp[2] = manl >> (GMP_LIMB_BITS - sc);
  178.       rp[1] = (manl << sc - GMP_NAIL_BITS) & GMP_NUMB_MASK;
  179.       if (sc >= 2 * GMP_NAIL_BITS)
  180. rp[0] = 0;
  181.       else
  182. rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
  183.     }
  184.   else
  185.     {
  186.       if (sc == 0)
  187. {
  188.   rp[2] = manl >> GMP_NAIL_BITS;
  189.   rp[1] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS) & GMP_NUMB_MASK;
  190.   rp[0] = 0;
  191.   exp--;
  192. }
  193.       else
  194. {
  195.   rp[2] = manl >> (GMP_LIMB_BITS - sc);
  196.   rp[1] = (manl >> GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
  197.   rp[0] = (manl << GMP_NUMB_BITS - GMP_NAIL_BITS + sc) & GMP_NUMB_MASK;
  198. }
  199.     }
  200. #endif
  201. #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE == 3
  202. #if GMP_NAIL_BITS == 0
  203.   if (sc != 0)
  204.     {
  205.       rp[2] = manh >> (GMP_LIMB_BITS - sc);
  206.       rp[1] = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
  207.       rp[0] = manl << sc;
  208.     }
  209.   else
  210.     {
  211.       rp[2] = manh;
  212.       rp[1] = manl;
  213.       rp[0] = 0;
  214.       exp--;
  215.     }
  216. #else
  217.   if (sc > GMP_NAIL_BITS)
  218.     {
  219.       rp[2] = (manh >> (GMP_LIMB_BITS - sc));
  220.       rp[1] = ((manh << (sc - GMP_NAIL_BITS)) |
  221.        (manl >> (GMP_LIMB_BITS - sc + GMP_NAIL_BITS))) & GMP_NUMB_MASK;
  222.       if (sc >= 2 * GMP_NAIL_BITS)
  223. rp[0] = (manl << sc - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
  224.       else
  225. rp[0] = manl >> (2 * GMP_NAIL_BITS - sc) & GMP_NUMB_MASK;
  226.     }
  227.   else
  228.     {
  229.       if (sc == 0)
  230. {
  231.   rp[2] = manh >> GMP_NAIL_BITS;
  232.   rp[1] = ((manh << GMP_NUMB_BITS - GMP_NAIL_BITS) | (manl >> 2 * GMP_NAIL_BITS)) & GMP_NUMB_MASK;
  233.   rp[0] = (manl << GMP_NUMB_BITS - 2 * GMP_NAIL_BITS) & GMP_NUMB_MASK;
  234.   exp--;
  235. }
  236.       else
  237. {
  238.   rp[2] = (manh >> (GMP_LIMB_BITS - sc));
  239.   rp[1] = (manh >> (GMP_NAIL_BITS - sc)) & GMP_NUMB_MASK;
  240.   rp[0] = ((manh << (GMP_NUMB_BITS - GMP_NAIL_BITS + sc))
  241.    | (manl >> (GMP_LIMB_BITS - (GMP_NUMB_BITS - GMP_NAIL_BITS + sc)))) & GMP_NUMB_MASK;
  242. }
  243.     }
  244. #endif
  245. #endif
  246. #if BITS_PER_PART == 32 && LIMBS_PER_DOUBLE > 3
  247.   if (sc == 0)
  248.     {
  249.       int i;
  250.       for (i = LIMBS_PER_DOUBLE - 1; i >= 0; i--)
  251. {
  252.   rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
  253.   manh = ((manh << GMP_NUMB_BITS)
  254.   | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
  255.   manl = manl << GMP_NUMB_BITS;
  256. }
  257.       exp--;
  258.     }
  259.   else
  260.     {
  261.       int i;
  262.       rp[LIMBS_PER_DOUBLE - 1] = (manh >> (GMP_LIMB_BITS - sc));
  263.       manh = (manh << sc) | (manl >> (GMP_LIMB_BITS - sc));
  264.       manl = (manl << sc);
  265.       for (i = LIMBS_PER_DOUBLE - 2; i >= 0; i--)
  266. {
  267.   rp[i] = manh >> (BITS_PER_ULONG - GMP_NUMB_BITS);
  268.   manh = ((manh << GMP_NUMB_BITS)
  269.   | (manl >> (BITS_PER_ULONG - GMP_NUMB_BITS)));
  270.   manl = manl << GMP_NUMB_BITS;
  271. }
  272.   }
  273. #endif
  274.   return exp;
  275. }