Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /****************************************************************
  2.  *
  3.  * The author of this software is David M. Gay.
  4.  *
  5.  * Copyright (c) 1991 by AT&T.
  6.  *
  7.  * Permission to use, copy, modify, and distribute this software for any
  8.  * purpose without fee is hereby granted, provided that this entire notice
  9.  * is included in all copies of any software which is or includes a copy
  10.  * or modification of this software and in all copies of the supporting
  11.  * documentation for such software.
  12.  *
  13.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17.  *
  18.  ***************************************************************/
  19.  
  20. /* Please send bug reports to
  21.         David M. Gay
  22.         AT&T Bell Laboratories, Room 2C-463
  23.         600 Mountain Avenue
  24.         Murray Hill, NJ 07974-2070
  25.         U.S.A.
  26.         dmg@research.att.com or research!dmg
  27.  */
  28.  
  29. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  30.  *
  31.  * This strtod returns a nearest machine number to the input decimal
  32.  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  33.  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  34.  * biased rounding (add half and chop).
  35.  *
  36.  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  37.  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  38.  *
  39.  * Modifications:
  40.  *
  41.  *      1. We only require IEEE, IBM, or VAX double-precision
  42.  *              arithmetic (not IEEE double-extended).
  43.  *      2. We get by with floating-point arithmetic in a case that
  44.  *              Clinger missed -- when we're computing d * 10^n
  45.  *              for a small integer d and the integer n is not too
  46.  *              much larger than 22 (the maximum integer k for which
  47.  *              we can represent 10^k exactly), we may be able to
  48.  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  49.  *      3. Rather than a bit-at-a-time adjustment of the binary
  50.  *              result in the hard case, we use floating-point
  51.  *              arithmetic to determine the adjustment to within
  52.  *              one bit; only in really hard cases do we need to
  53.  *              compute a second residual.
  54.  *      4. Because of 3., we don't need a large table of powers of 10
  55.  *              for ten-to-e (just some small tables, e.g. of 10^k
  56.  *              for 0 <= k <= 22).
  57.  */
  58.  
  59. /*
  60.  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  61.  *      significant byte has the lowest address.
  62.  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  63.  *      significant byte has the lowest address.
  64.  * #define Sudden_Underflow for IEEE-format machines without gradual
  65.  *      underflow (i.e., that flush to zero on underflow).
  66.  * #define IBM for IBM mainframe-style floating-point arithmetic.
  67.  * #define VAX for VAX-style floating-point arithmetic.
  68.  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
  69.  * #define No_leftright to omit left-right logic in fast floating-point
  70.  *      computation of dtoa.
  71.  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  72.  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  73.  *      that use extended-precision instructions to compute rounded
  74.  *      products and quotients) with IBM.
  75.  * #define ROUND_BIASED for IEEE-format with biased rounding.
  76.  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  77.  *      products but inaccurate quotients, e.g., for Intel i860.
  78.  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
  79.  *      integer arithmetic.  Whether this speeds things up or slows things
  80.  *      down depends on the machine and the number being converted.
  81.  */
  82.  
  83. #include <_ansi.h>
  84. #include <stdlib.h>
  85. #include <string.h>
  86. #include <reent.h>
  87. #include "mprec.h"
  88.  
  89. /* This is defined in sys/reent.h as (sizeof (size_t) << 3) now, as in NetBSD.
  90.    The old value of 15 was wrong and made newlib vulnerable against buffer
  91.    overrun attacks (CVE-2009-0689), same as other implementations of gdtoa
  92.    based on BSD code.
  93. #define _Kmax 15
  94. */
  95.  
  96. _Bigint *
  97. _DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k)
  98. {
  99.   int x;
  100.   _Bigint *rv ;
  101.  
  102.   _REENT_CHECK_MP(ptr);
  103.   if (_REENT_MP_FREELIST(ptr) == NULL)
  104.     {
  105.       /* Allocate a list of pointers to the mprec objects */
  106.       _REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr,
  107.                                                       sizeof (struct _Bigint *),
  108.                                                       _Kmax + 1);
  109.       if (_REENT_MP_FREELIST(ptr) == NULL)
  110.         {
  111.           return NULL;
  112.         }
  113.     }
  114.  
  115.   if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0)
  116.     {
  117.       _REENT_MP_FREELIST(ptr)[k] = rv->_next;
  118.     }
  119.   else
  120.     {
  121.       x = 1 << k;
  122.       /* Allocate an mprec Bigint and stick in in the freelist */
  123.       rv = (_Bigint *) _calloc_r (ptr,
  124.                                   1,
  125.                                   sizeof (_Bigint) +
  126.                                   (x-1) * sizeof(rv->_x));
  127.       if (rv == NULL) return NULL;
  128.       rv->_k = k;
  129.       rv->_maxwds = x;
  130.     }
  131.   rv->_sign = rv->_wds = 0;
  132.   return rv;
  133. }
  134.  
  135. void
  136. _DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v)
  137. {
  138.   _REENT_CHECK_MP(ptr);
  139.   if (v)
  140.     {
  141.       v->_next = _REENT_MP_FREELIST(ptr)[v->_k];
  142.       _REENT_MP_FREELIST(ptr)[v->_k] = v;
  143.     }
  144. }
  145.  
  146. _Bigint *
  147. _DEFUN (multadd, (ptr, b, m, a),
  148.         struct _reent *ptr _AND
  149.         _Bigint * b _AND
  150.         int m _AND
  151.         int a)
  152. {
  153.   int i, wds;
  154.   __ULong *x, y;
  155. #ifdef Pack_32
  156.   __ULong xi, z;
  157. #endif
  158.   _Bigint *b1;
  159.  
  160.   wds = b->_wds;
  161.   x = b->_x;
  162.   i = 0;
  163.   do
  164.     {
  165. #ifdef Pack_32
  166.       xi = *x;
  167.       y = (xi & 0xffff) * m + a;
  168.       z = (xi >> 16) * m + (y >> 16);
  169.       a = (int) (z >> 16);
  170.       *x++ = (z << 16) + (y & 0xffff);
  171. #else
  172.       y = *x * m + a;
  173.       a = (int) (y >> 16);
  174.       *x++ = y & 0xffff;
  175. #endif
  176.     }
  177.   while (++i < wds);
  178.   if (a)
  179.     {
  180.       if (wds >= b->_maxwds)
  181.         {
  182.           b1 = Balloc (ptr, b->_k + 1);
  183.           Bcopy (b1, b);
  184.           Bfree (ptr, b);
  185.           b = b1;
  186.         }
  187.       b->_x[wds++] = a;
  188.       b->_wds = wds;
  189.     }
  190.   return b;
  191. }
  192.  
  193. _Bigint *
  194. _DEFUN (s2b, (ptr, s, nd0, nd, y9),
  195.         struct _reent * ptr _AND
  196.         _CONST char *s _AND
  197.         int nd0 _AND
  198.         int nd _AND
  199.         __ULong y9)
  200. {
  201.   _Bigint *b;
  202.   int i, k;
  203.   __Long x, y;
  204.  
  205.   x = (nd + 8) / 9;
  206.   for (k = 0, y = 1; x > y; y <<= 1, k++);
  207. #ifdef Pack_32
  208.   b = Balloc (ptr, k);
  209.   b->_x[0] = y9;
  210.   b->_wds = 1;
  211. #else
  212.   b = Balloc (ptr, k + 1);
  213.   b->_x[0] = y9 & 0xffff;
  214.   b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
  215. #endif
  216.  
  217.   i = 9;
  218.   if (9 < nd0)
  219.     {
  220.       s += 9;
  221.       do
  222.         b = multadd (ptr, b, 10, *s++ - '0');
  223.       while (++i < nd0);
  224.       s++;
  225.     }
  226.   else
  227.     s += 10;
  228.   for (; i < nd; i++)
  229.     b = multadd (ptr, b, 10, *s++ - '0');
  230.   return b;
  231. }
  232.  
  233. int
  234. _DEFUN (hi0bits,
  235.         (x), register __ULong x)
  236. {
  237.   register int k = 0;
  238.  
  239.   if (!(x & 0xffff0000))
  240.     {
  241.       k = 16;
  242.       x <<= 16;
  243.     }
  244.   if (!(x & 0xff000000))
  245.     {
  246.       k += 8;
  247.       x <<= 8;
  248.     }
  249.   if (!(x & 0xf0000000))
  250.     {
  251.       k += 4;
  252.       x <<= 4;
  253.     }
  254.   if (!(x & 0xc0000000))
  255.     {
  256.       k += 2;
  257.       x <<= 2;
  258.     }
  259.   if (!(x & 0x80000000))
  260.     {
  261.       k++;
  262.       if (!(x & 0x40000000))
  263.         return 32;
  264.     }
  265.   return k;
  266. }
  267.  
  268. int
  269. _DEFUN (lo0bits, (y), __ULong *y)
  270. {
  271.   register int k;
  272.   register __ULong x = *y;
  273.  
  274.   if (x & 7)
  275.     {
  276.       if (x & 1)
  277.         return 0;
  278.       if (x & 2)
  279.         {
  280.           *y = x >> 1;
  281.           return 1;
  282.         }
  283.       *y = x >> 2;
  284.       return 2;
  285.     }
  286.   k = 0;
  287.   if (!(x & 0xffff))
  288.     {
  289.       k = 16;
  290.       x >>= 16;
  291.     }
  292.   if (!(x & 0xff))
  293.     {
  294.       k += 8;
  295.       x >>= 8;
  296.     }
  297.   if (!(x & 0xf))
  298.     {
  299.       k += 4;
  300.       x >>= 4;
  301.     }
  302.   if (!(x & 0x3))
  303.     {
  304.       k += 2;
  305.       x >>= 2;
  306.     }
  307.   if (!(x & 1))
  308.     {
  309.       k++;
  310.       x >>= 1;
  311.       if (!x & 1)
  312.         return 32;
  313.     }
  314.   *y = x;
  315.   return k;
  316. }
  317.  
  318. _Bigint *
  319. _DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i)
  320. {
  321.   _Bigint *b;
  322.  
  323.   b = Balloc (ptr, 1);
  324.   b->_x[0] = i;
  325.   b->_wds = 1;
  326.   return b;
  327. }
  328.  
  329. _Bigint *
  330. _DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b)
  331. {
  332.   _Bigint *c;
  333.   int k, wa, wb, wc;
  334.   __ULong carry, y, z;
  335.   __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  336. #ifdef Pack_32
  337.   __ULong z2;
  338. #endif
  339.  
  340.   if (a->_wds < b->_wds)
  341.     {
  342.       c = a;
  343.       a = b;
  344.       b = c;
  345.     }
  346.   k = a->_k;
  347.   wa = a->_wds;
  348.   wb = b->_wds;
  349.   wc = wa + wb;
  350.   if (wc > a->_maxwds)
  351.     k++;
  352.   c = Balloc (ptr, k);
  353.   for (x = c->_x, xa = x + wc; x < xa; x++)
  354.     *x = 0;
  355.   xa = a->_x;
  356.   xae = xa + wa;
  357.   xb = b->_x;
  358.   xbe = xb + wb;
  359.   xc0 = c->_x;
  360. #ifdef Pack_32
  361.   for (; xb < xbe; xb++, xc0++)
  362.     {
  363.       if ((y = *xb & 0xffff) != 0)
  364.         {
  365.           x = xa;
  366.           xc = xc0;
  367.           carry = 0;
  368.           do
  369.             {
  370.               z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  371.               carry = z >> 16;
  372.               z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  373.               carry = z2 >> 16;
  374.               Storeinc (xc, z2, z);
  375.             }
  376.           while (x < xae);
  377.           *xc = carry;
  378.         }
  379.       if ((y = *xb >> 16) != 0)
  380.         {
  381.           x = xa;
  382.           xc = xc0;
  383.           carry = 0;
  384.           z2 = *xc;
  385.           do
  386.             {
  387.               z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  388.               carry = z >> 16;
  389.               Storeinc (xc, z, z2);
  390.               z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  391.               carry = z2 >> 16;
  392.             }
  393.           while (x < xae);
  394.           *xc = z2;
  395.         }
  396.     }
  397. #else
  398.   for (; xb < xbe; xc0++)
  399.     {
  400.       if (y = *xb++)
  401.         {
  402.           x = xa;
  403.           xc = xc0;
  404.           carry = 0;
  405.           do
  406.             {
  407.               z = *x++ * y + *xc + carry;
  408.               carry = z >> 16;
  409.               *xc++ = z & 0xffff;
  410.             }
  411.           while (x < xae);
  412.           *xc = carry;
  413.         }
  414.     }
  415. #endif
  416.   for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
  417.   c->_wds = wc;
  418.   return c;
  419. }
  420.  
  421. _Bigint *
  422. _DEFUN (pow5mult,
  423.         (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
  424. {
  425.   _Bigint *b1, *p5, *p51;
  426.   int i;
  427.   static _CONST int p05[3] = {5, 25, 125};
  428.  
  429.   if ((i = k & 3) != 0)
  430.     b = multadd (ptr, b, p05[i - 1], 0);
  431.  
  432.   if (!(k >>= 2))
  433.     return b;
  434.   _REENT_CHECK_MP(ptr);
  435.   if (!(p5 = _REENT_MP_P5S(ptr)))
  436.     {
  437.       /* first time */
  438.       p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625);
  439.       p5->_next = 0;
  440.     }
  441.   for (;;)
  442.     {
  443.       if (k & 1)
  444.         {
  445.           b1 = mult (ptr, b, p5);
  446.           Bfree (ptr, b);
  447.           b = b1;
  448.         }
  449.       if (!(k >>= 1))
  450.         break;
  451.       if (!(p51 = p5->_next))
  452.         {
  453.           p51 = p5->_next = mult (ptr, p5, p5);
  454.           p51->_next = 0;
  455.         }
  456.       p5 = p51;
  457.     }
  458.   return b;
  459. }
  460.  
  461. _Bigint *
  462. _DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
  463. {
  464.   int i, k1, n, n1;
  465.   _Bigint *b1;
  466.   __ULong *x, *x1, *xe, z;
  467.  
  468. #ifdef Pack_32
  469.   n = k >> 5;
  470. #else
  471.   n = k >> 4;
  472. #endif
  473.   k1 = b->_k;
  474.   n1 = n + b->_wds + 1;
  475.   for (i = b->_maxwds; n1 > i; i <<= 1)
  476.     k1++;
  477.   b1 = Balloc (ptr, k1);
  478.   x1 = b1->_x;
  479.   for (i = 0; i < n; i++)
  480.     *x1++ = 0;
  481.   x = b->_x;
  482.   xe = x + b->_wds;
  483. #ifdef Pack_32
  484.   if (k &= 0x1f)
  485.     {
  486.       k1 = 32 - k;
  487.       z = 0;
  488.       do
  489.         {
  490.           *x1++ = *x << k | z;
  491.           z = *x++ >> k1;
  492.         }
  493.       while (x < xe);
  494.       if ((*x1 = z) != 0)
  495.         ++n1;
  496.     }
  497. #else
  498.   if (k &= 0xf)
  499.     {
  500.       k1 = 16 - k;
  501.       z = 0;
  502.       do
  503.         {
  504.           *x1++ = *x << k & 0xffff | z;
  505.           z = *x++ >> k1;
  506.         }
  507.       while (x < xe);
  508.       if (*x1 = z)
  509.         ++n1;
  510.     }
  511. #endif
  512.   else
  513.     do
  514.       *x1++ = *x++;
  515.     while (x < xe);
  516.   b1->_wds = n1 - 1;
  517.   Bfree (ptr, b);
  518.   return b1;
  519. }
  520.  
  521. int
  522. _DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b)
  523. {
  524.   __ULong *xa, *xa0, *xb, *xb0;
  525.   int i, j;
  526.  
  527.   i = a->_wds;
  528.   j = b->_wds;
  529. #ifdef DEBUG
  530.   if (i > 1 && !a->_x[i - 1])
  531.     Bug ("cmp called with a->_x[a->_wds-1] == 0");
  532.   if (j > 1 && !b->_x[j - 1])
  533.     Bug ("cmp called with b->_x[b->_wds-1] == 0");
  534. #endif
  535.   if (i -= j)
  536.     return i;
  537.   xa0 = a->_x;
  538.   xa = xa0 + j;
  539.   xb0 = b->_x;
  540.   xb = xb0 + j;
  541.   for (;;)
  542.     {
  543.       if (*--xa != *--xb)
  544.         return *xa < *xb ? -1 : 1;
  545.       if (xa <= xa0)
  546.         break;
  547.     }
  548.   return 0;
  549. }
  550.  
  551. _Bigint *
  552. _DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND
  553.         _Bigint * a _AND _Bigint * b)
  554. {
  555.   _Bigint *c;
  556.   int i, wa, wb;
  557.   __Long borrow, y;             /* We need signed shifts here. */
  558.   __ULong *xa, *xae, *xb, *xbe, *xc;
  559. #ifdef Pack_32
  560.   __Long z;
  561. #endif
  562.  
  563.   i = cmp (a, b);
  564.   if (!i)
  565.     {
  566.       c = Balloc (ptr, 0);
  567.       c->_wds = 1;
  568.       c->_x[0] = 0;
  569.       return c;
  570.     }
  571.   if (i < 0)
  572.     {
  573.       c = a;
  574.       a = b;
  575.       b = c;
  576.       i = 1;
  577.     }
  578.   else
  579.     i = 0;
  580.   c = Balloc (ptr, a->_k);
  581.   c->_sign = i;
  582.   wa = a->_wds;
  583.   xa = a->_x;
  584.   xae = xa + wa;
  585.   wb = b->_wds;
  586.   xb = b->_x;
  587.   xbe = xb + wb;
  588.   xc = c->_x;
  589.   borrow = 0;
  590. #ifdef Pack_32
  591.   do
  592.     {
  593.       y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
  594.       borrow = y >> 16;
  595.       Sign_Extend (borrow, y);
  596.       z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
  597.       borrow = z >> 16;
  598.       Sign_Extend (borrow, z);
  599.       Storeinc (xc, z, y);
  600.     }
  601.   while (xb < xbe);
  602.   while (xa < xae)
  603.     {
  604.       y = (*xa & 0xffff) + borrow;
  605.       borrow = y >> 16;
  606.       Sign_Extend (borrow, y);
  607.       z = (*xa++ >> 16) + borrow;
  608.       borrow = z >> 16;
  609.       Sign_Extend (borrow, z);
  610.       Storeinc (xc, z, y);
  611.     }
  612. #else
  613.   do
  614.     {
  615.       y = *xa++ - *xb++ + borrow;
  616.       borrow = y >> 16;
  617.       Sign_Extend (borrow, y);
  618.       *xc++ = y & 0xffff;
  619.     }
  620.   while (xb < xbe);
  621.   while (xa < xae)
  622.     {
  623.       y = *xa++ + borrow;
  624.       borrow = y >> 16;
  625.       Sign_Extend (borrow, y);
  626.       *xc++ = y & 0xffff;
  627.     }
  628. #endif
  629.   while (!*--xc)
  630.     wa--;
  631.   c->_wds = wa;
  632.   return c;
  633. }
  634.  
  635. double
  636. _DEFUN (ulp, (_x), double _x)
  637. {
  638.   union double_union x, a;
  639.   register __Long L;
  640.  
  641.   x.d = _x;
  642.  
  643.   L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
  644. #ifndef Sudden_Underflow
  645.   if (L > 0)
  646.     {
  647. #endif
  648. #ifdef IBM
  649.       L |= Exp_msk1 >> 4;
  650. #endif
  651.       word0 (a) = L;
  652. #ifndef _DOUBLE_IS_32BITS
  653.       word1 (a) = 0;
  654. #endif
  655.  
  656. #ifndef Sudden_Underflow
  657.     }
  658.   else
  659.     {
  660.       L = -L >> Exp_shift;
  661.       if (L < Exp_shift)
  662.         {
  663.           word0 (a) = 0x80000 >> L;
  664. #ifndef _DOUBLE_IS_32BITS
  665.           word1 (a) = 0;
  666. #endif
  667.         }
  668.       else
  669.         {
  670.           word0 (a) = 0;
  671.           L -= Exp_shift;
  672. #ifndef _DOUBLE_IS_32BITS
  673.          word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
  674. #endif
  675.         }
  676.     }
  677. #endif
  678.   return a.d;
  679. }
  680.  
  681. double
  682. _DEFUN (b2d, (a, e),
  683.         _Bigint * a _AND int *e)
  684. {
  685.   __ULong *xa, *xa0, w, y, z;
  686.   int k;
  687.   union double_union d;
  688. #ifdef VAX
  689.   __ULong d0, d1;
  690. #else
  691. #define d0 word0(d)
  692. #define d1 word1(d)
  693. #endif
  694.  
  695.   xa0 = a->_x;
  696.   xa = xa0 + a->_wds;
  697.   y = *--xa;
  698. #ifdef DEBUG
  699.   if (!y)
  700.     Bug ("zero y in b2d");
  701. #endif
  702.   k = hi0bits (y);
  703.   *e = 32 - k;
  704. #ifdef Pack_32
  705.   if (k < Ebits)
  706.     {
  707.       d0 = Exp_1 | y >> (Ebits - k);
  708.       w = xa > xa0 ? *--xa : 0;
  709. #ifndef _DOUBLE_IS_32BITS
  710.       d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
  711. #endif
  712.       goto ret_d;
  713.     }
  714.   z = xa > xa0 ? *--xa : 0;
  715.   if (k -= Ebits)
  716.     {
  717.       d0 = Exp_1 | y << k | z >> (32 - k);
  718.       y = xa > xa0 ? *--xa : 0;
  719. #ifndef _DOUBLE_IS_32BITS
  720.       d1 = z << k | y >> (32 - k);
  721. #endif
  722.     }
  723.   else
  724.     {
  725.       d0 = Exp_1 | y;
  726. #ifndef _DOUBLE_IS_32BITS
  727.       d1 = z;
  728. #endif
  729.     }
  730. #else
  731.   if (k < Ebits + 16)
  732.     {
  733.       z = xa > xa0 ? *--xa : 0;
  734.       d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  735.       w = xa > xa0 ? *--xa : 0;
  736.       y = xa > xa0 ? *--xa : 0;
  737.       d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  738.       goto ret_d;
  739.     }
  740.   z = xa > xa0 ? *--xa : 0;
  741.   w = xa > xa0 ? *--xa : 0;
  742.   k -= Ebits + 16;
  743.   d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  744.   y = xa > xa0 ? *--xa : 0;
  745.   d1 = w << k + 16 | y << k;
  746. #endif
  747. ret_d:
  748. #ifdef VAX
  749.   word0 (d) = d0 >> 16 | d0 << 16;
  750.   word1 (d) = d1 >> 16 | d1 << 16;
  751. #else
  752. #undef d0
  753. #undef d1
  754. #endif
  755.   return d.d;
  756. }
  757.  
  758. _Bigint *
  759. _DEFUN (d2b,
  760.         (ptr, _d, e, bits),
  761.         struct _reent * ptr _AND
  762.         double _d _AND
  763.         int *e _AND
  764.         int *bits)
  765.  
  766. {
  767.   union double_union d;
  768.   _Bigint *b;
  769.   int de, i, k;
  770.   __ULong *x, y, z;
  771. #ifdef VAX
  772.   __ULong d0, d1;
  773. #endif
  774.   d.d = _d;
  775. #ifdef VAX
  776.   d0 = word0 (d) >> 16 | word0 (d) << 16;
  777.   d1 = word1 (d) >> 16 | word1 (d) << 16;
  778. #else
  779. #define d0 word0(d)
  780. #define d1 word1(d)
  781.   d.d = _d;
  782. #endif
  783.  
  784. #ifdef Pack_32
  785.   b = Balloc (ptr, 1);
  786. #else
  787.   b = Balloc (ptr, 2);
  788. #endif
  789.   x = b->_x;
  790.  
  791.   z = d0 & Frac_mask;
  792.   d0 &= 0x7fffffff;             /* clear sign bit, which we ignore */
  793. #ifdef Sudden_Underflow
  794.   de = (int) (d0 >> Exp_shift);
  795. #ifndef IBM
  796.   z |= Exp_msk11;
  797. #endif
  798. #else
  799.   if ((de = (int) (d0 >> Exp_shift)) != 0)
  800.     z |= Exp_msk1;
  801. #endif
  802. #ifdef Pack_32
  803. #ifndef _DOUBLE_IS_32BITS
  804.   if (d1)
  805.     {
  806.       y = d1;
  807.       k = lo0bits (&y);
  808.       if (k)
  809.         {
  810.          x[0] = y | z << (32 - k);
  811.           z >>= k;
  812.         }
  813.       else
  814.         x[0] = y;
  815.       i = b->_wds = (x[1] = z) ? 2 : 1;
  816.     }
  817.   else
  818. #endif
  819.     {
  820. #ifdef DEBUG
  821.       if (!z)
  822.         Bug ("Zero passed to d2b");
  823. #endif
  824.       k = lo0bits (&z);
  825.       x[0] = z;
  826.       i = b->_wds = 1;
  827. #ifndef _DOUBLE_IS_32BITS
  828.       k += 32;
  829. #endif
  830.     }
  831. #else
  832.   if (d1)
  833.     {
  834.       y = d1;
  835.       k = lo0bits (&y);
  836.       if (k)
  837.         if (k >= 16)
  838.           {
  839.             x[0] = y | z << 32 - k & 0xffff;
  840.             x[1] = z >> k - 16 & 0xffff;
  841.             x[2] = z >> k;
  842.             i = 2;
  843.           }
  844.         else
  845.           {
  846.             x[0] = y & 0xffff;
  847.             x[1] = y >> 16 | z << 16 - k & 0xffff;
  848.             x[2] = z >> k & 0xffff;
  849.             x[3] = z >> k + 16;
  850.             i = 3;
  851.           }
  852.       else
  853.         {
  854.           x[0] = y & 0xffff;
  855.           x[1] = y >> 16;
  856.           x[2] = z & 0xffff;
  857.           x[3] = z >> 16;
  858.           i = 3;
  859.         }
  860.     }
  861.   else
  862.     {
  863. #ifdef DEBUG
  864.       if (!z)
  865.         Bug ("Zero passed to d2b");
  866. #endif
  867.       k = lo0bits (&z);
  868.       if (k >= 16)
  869.         {
  870.           x[0] = z;
  871.           i = 0;
  872.         }
  873.       else
  874.         {
  875.           x[0] = z & 0xffff;
  876.           x[1] = z >> 16;
  877.           i = 1;
  878.         }
  879.       k += 32;
  880.     }
  881.   while (!x[i])
  882.     --i;
  883.   b->_wds = i + 1;
  884. #endif
  885. #ifndef Sudden_Underflow
  886.   if (de)
  887.     {
  888. #endif
  889. #ifdef IBM
  890.       *e = (de - Bias - (P - 1) << 2) + k;
  891.       *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
  892. #else
  893.       *e = de - Bias - (P - 1) + k;
  894.       *bits = P - k;
  895. #endif
  896. #ifndef Sudden_Underflow
  897.     }
  898.   else
  899.     {
  900.       *e = de - Bias - (P - 1) + 1 + k;
  901. #ifdef Pack_32
  902.       *bits = 32 * i - hi0bits (x[i - 1]);
  903. #else
  904.       *bits = (i + 2) * 16 - hi0bits (x[i]);
  905. #endif
  906.     }
  907. #endif
  908.   return b;
  909. }
  910. #undef d0
  911. #undef d1
  912.  
  913. double
  914. _DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b)
  915.  
  916. {
  917.   union double_union da, db;
  918.   int k, ka, kb;
  919.  
  920.   da.d = b2d (a, &ka);
  921.   db.d = b2d (b, &kb);
  922. #ifdef Pack_32
  923.   k = ka - kb + 32 * (a->_wds - b->_wds);
  924. #else
  925.   k = ka - kb + 16 * (a->_wds - b->_wds);
  926. #endif
  927. #ifdef IBM
  928.   if (k > 0)
  929.     {
  930.       word0 (da) += (k >> 2) * Exp_msk1;
  931.       if (k &= 3)
  932.         da.d *= 1 << k;
  933.     }
  934.   else
  935.     {
  936.       k = -k;
  937.       word0 (db) += (k >> 2) * Exp_msk1;
  938.       if (k &= 3)
  939.         db.d *= 1 << k;
  940.     }
  941. #else
  942.   if (k > 0)
  943.     word0 (da) += k * Exp_msk1;
  944.   else
  945.     {
  946.       k = -k;
  947.       word0 (db) += k * Exp_msk1;
  948.     }
  949. #endif
  950.   return da.d / db.d;
  951. }
  952.  
  953.  
  954. _CONST double
  955.   tens[] =
  956. {
  957.   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  958.   1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  959.   1e20, 1e21, 1e22, 1e23, 1e24
  960.  
  961. };
  962.  
  963. #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
  964. _CONST double bigtens[] =
  965. {1e16, 1e32, 1e64, 1e128, 1e256};
  966.  
  967. _CONST double tinytens[] =
  968. {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
  969. #else
  970. _CONST double bigtens[] =
  971. {1e16, 1e32};
  972.  
  973. _CONST double tinytens[] =
  974. {1e-16, 1e-32};
  975. #endif
  976.  
  977.  
  978. double
  979. _DEFUN (_mprec_log10, (dig),
  980.         int dig)
  981. {
  982.   double v = 1.0;
  983.   if (dig < 24)
  984.     return tens[dig];
  985.   while (dig > 0)
  986.     {
  987.       v *= 10;
  988.       dig--;
  989.     }
  990.   return v;
  991. }
  992.  
  993. void
  994. _DEFUN (copybits, (c, n, b),
  995.         __ULong *c _AND
  996.         int n _AND
  997.         _Bigint *b)
  998. {
  999.         __ULong *ce, *x, *xe;
  1000. #ifdef Pack_16
  1001.         int nw, nw1;
  1002. #endif
  1003.  
  1004.         ce = c + ((n-1) >> kshift) + 1;
  1005.         x = b->_x;
  1006. #ifdef Pack_32
  1007.         xe = x + b->_wds;
  1008.         while(x < xe)
  1009.                 *c++ = *x++;
  1010. #else
  1011.         nw = b->_wds;
  1012.         nw1 = nw & 1;
  1013.         for(xe = x + (nw - nw1); x < xe; x += 2)
  1014.                 Storeinc(c, x[1], x[0]);
  1015.         if (nw1)
  1016.                 *c++ = *x;
  1017. #endif
  1018.         while(c < ce)
  1019.                 *c++ = 0;
  1020. }
  1021.  
  1022. __ULong
  1023. _DEFUN (any_on, (b, k),
  1024.         _Bigint *b _AND
  1025.         int k)
  1026. {
  1027.         int n, nwds;
  1028.         __ULong *x, *x0, x1, x2;
  1029.  
  1030.         x = b->_x;
  1031.         nwds = b->_wds;
  1032.         n = k >> kshift;
  1033.         if (n > nwds)
  1034.                 n = nwds;
  1035.         else if (n < nwds && (k &= kmask)) {
  1036.                 x1 = x2 = x[n];
  1037.                 x1 >>= k;
  1038.                 x1 <<= k;
  1039.                 if (x1 != x2)
  1040.                         return 1;
  1041.                 }
  1042.         x0 = x;
  1043.         x += n;
  1044.         while(x > x0)
  1045.                 if (*--x)
  1046.                         return 1;
  1047.         return 0;
  1048. }
  1049.  
  1050.