Contiki 2.6

small_mprec.c

00001 /****************************************************************
00002  *
00003  * The author of this software is David M. Gay.
00004  *
00005  * Copyright (c) 1991 by AT&T.
00006  *
00007  * Permission to use, copy, modify, and distribute this software for any
00008  * purpose without fee is hereby granted, provided that this entire notice
00009  * is included in all copies of any software which is or includes a copy
00010  * or modification of this software and in all copies of the supporting
00011  * documentation for such software.
00012  *
00013  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00014  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
00015  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00016  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00017  *
00018  ***************************************************************/
00019 
00020 /* Please send bug reports to
00021         David M. Gay
00022         AT&T Bell Laboratories, Room 2C-463
00023         600 Mountain Avenue
00024         Murray Hill, NJ 07974-2070
00025         U.S.A.
00026         dmg@research.att.com or research!dmg
00027  */
00028 
00029 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
00030  *
00031  * This strtod returns a nearest machine number to the input decimal
00032  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
00033  * broken by the IEEE round-even rule.  Otherwise ties are broken by
00034  * biased rounding (add half and chop).
00035  *
00036  * Inspired loosely by William D. Clinger's paper "How to Read Floating
00037  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
00038  *
00039  * Modifications:
00040  *
00041  *      1. We only require IEEE, IBM, or VAX double-precision
00042  *              arithmetic (not IEEE double-extended).
00043  *      2. We get by with floating-point arithmetic in a case that
00044  *              Clinger missed -- when we're computing d * 10^n
00045  *              for a small integer d and the integer n is not too
00046  *              much larger than 22 (the maximum integer k for which
00047  *              we can represent 10^k exactly), we may be able to
00048  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
00049  *      3. Rather than a bit-at-a-time adjustment of the binary
00050  *              result in the hard case, we use floating-point
00051  *              arithmetic to determine the adjustment to within
00052  *              one bit; only in really hard cases do we need to
00053  *              compute a second residual.
00054  *      4. Because of 3., we don't need a large table of powers of 10
00055  *              for ten-to-e (just some small tables, e.g. of 10^k
00056  *              for 0 <= k <= 22).
00057  */
00058 
00059 /*
00060  * #define IEEE_8087 for IEEE-arithmetic machines where the least
00061  *      significant byte has the lowest address.
00062  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
00063  *      significant byte has the lowest address.
00064  * #define Sudden_Underflow for IEEE-format machines without gradual
00065  *      underflow (i.e., that flush to zero on underflow).
00066  * #define IBM for IBM mainframe-style floating-point arithmetic.
00067  * #define VAX for VAX-style floating-point arithmetic.
00068  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
00069  * #define No_leftright to omit left-right logic in fast floating-point
00070  *      computation of dtoa.
00071  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
00072  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
00073  *      that use extended-precision instructions to compute rounded
00074  *      products and quotients) with IBM.
00075  * #define ROUND_BIASED for IEEE-format with biased rounding.
00076  * #define Inaccurate_Divide for IEEE-format with correctly rounded
00077  *      products but inaccurate quotients, e.g., for Intel i860.
00078  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
00079  *      integer arithmetic.  Whether this speeds things up or slows things
00080  *      down depends on the machine and the number being converted.
00081  */
00082 
00083 
00084 
00085 #include <_ansi.h>
00086 #include <stdlib.h>
00087 #include <string.h>
00088 #include <reent.h>
00089 #include "small_mprec.h"
00090 
00091 /* reent.c knows this value */
00092 #define _Kmax 15
00093 
00094 #if defined (_SMALL_PRINTF) || defined(SMALL_SCANF)
00095 #define SMALL_LIB
00096 #else
00097 #define FULL_LIB
00098 #endif
00099 
00100 /* SMALL_LIB is only defined if _SMALL_PRINTF or SMALL_SCANF have been defined
00101  * which means that if you do only specified _SMALL_PRINTF and not SMALL_SCANF
00102  * optimization about call of balloc in small_dtoa_r.c or smal_strtod.c will be
00103  * performed for both printf and scanf 
00104  */
00105 
00106 
00107 
00108 #ifdef FULL_LIB
00109 
00110 
00111 
00112 _Bigint *
00113 _DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k)
00114 {
00115   int x;
00116   _Bigint *rv ;
00117 
00118   _REENT_CHECK_MP(ptr);
00119   if (_REENT_MP_FREELIST(ptr) == NULL)
00120     {
00121       /* Allocate a list of pointers to the mprec objects */
00122       _REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr, 
00123                                                       sizeof (struct _Bigint *),
00124                                                       _Kmax + 1);
00125       if (_REENT_MP_FREELIST(ptr) == NULL)
00126         {
00127           return NULL;
00128         }
00129     }
00130 
00131   if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0)
00132     {
00133       _REENT_MP_FREELIST(ptr)[k] = rv->_next;
00134     }
00135   else
00136     {
00137       x = 1 << k;
00138       /* Allocate an mprec Bigint and stick in in the freelist */
00139       rv = (_Bigint *) _calloc_r (ptr,
00140                                   1,
00141                                   sizeof (_Bigint) +
00142                                   (x-1) * sizeof(rv->_x));
00143       if (rv == NULL) return NULL;
00144       rv->_k = k;
00145       rv->_maxwds = x;
00146     }
00147   rv->_sign = rv->_wds = 0;
00148   return rv;
00149 }
00150 
00151 void
00152 _DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v)
00153 {
00154   _REENT_CHECK_MP(ptr);
00155   if (v)
00156     {
00157       v->_next = _REENT_MP_FREELIST(ptr)[v->_k];
00158       _REENT_MP_FREELIST(ptr)[v->_k] = v;
00159     }
00160 }
00161 
00162 
00163 _Bigint *
00164 _DEFUN (multadd, (ptr, b, m, a),
00165         struct _reent *ptr _AND
00166         _Bigint * b _AND
00167         int m _AND
00168         int a)
00169 {
00170   int i, wds;
00171   __ULong *x, y;
00172 #ifdef Pack_32
00173   __ULong xi, z;
00174 #endif
00175   _Bigint *b1;
00176 
00177   wds = b->_wds;
00178   x = b->_x;
00179   i = 0;
00180   do
00181     {
00182 #ifdef Pack_32
00183       xi = *x;
00184       y = (xi & 0xffff) * m + a;
00185       z = (xi >> 16) * m + (y >> 16);
00186       a = (int) (z >> 16);
00187       *x++ = (z << 16) + (y & 0xffff);
00188 #else
00189       y = *x * m + a;
00190       a = (int) (y >> 16);
00191       *x++ = y & 0xffff;
00192 #endif
00193     }
00194   while (++i < wds);
00195   if (a)
00196     {
00197       if (wds >= b->_maxwds)
00198         {
00199           b1 = Balloc (ptr, b->_k + 1);
00200           Bcopy (b1, b);
00201           Bfree (ptr, b);
00202           b = b1;
00203         }
00204       b->_x[wds++] = a;
00205       b->_wds = wds;
00206     }
00207   return b;
00208 }
00209 
00210 _Bigint *
00211 _DEFUN (s2b, (ptr, s, nd0, nd, y9),
00212         struct _reent * ptr _AND
00213         _CONST char *s _AND
00214         int nd0 _AND
00215         int nd _AND
00216         __ULong y9)
00217 {
00218   _Bigint *b;
00219   int i, k;
00220   __Long x, y;
00221 
00222   x = (nd + 8) / 9;
00223   for (k = 0, y = 1; x > y; y <<= 1, k++);
00224 #ifdef Pack_32
00225   b = Balloc (ptr, k);
00226   b->_x[0] = y9;
00227   b->_wds = 1;
00228 #else
00229   b = Balloc (ptr, k + 1);
00230   b->_x[0] = y9 & 0xffff;
00231   b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
00232 #endif
00233 
00234   i = 9;
00235   if (9 < nd0)
00236     {
00237       s += 9;
00238       do
00239         b = multadd (ptr, b, 10, *s++ - '0');
00240       while (++i < nd0);
00241       s++;
00242     }
00243   else
00244     s += 10;
00245   for (; i < nd; i++)
00246     b = multadd (ptr, b, 10, *s++ - '0');
00247   return b;
00248 }
00249 
00250 int
00251 _DEFUN (hi0bits,
00252         (x), register __ULong x)
00253 {
00254   register int k = 0;
00255 
00256   if (!(x & 0xffff0000))
00257     {
00258       k = 16;
00259       x <<= 16;
00260     }
00261   if (!(x & 0xff000000))
00262     {
00263       k += 8;
00264       x <<= 8;
00265     }
00266   if (!(x & 0xf0000000))
00267     {
00268       k += 4;
00269       x <<= 4;
00270     }
00271   if (!(x & 0xc0000000))
00272     {
00273       k += 2;
00274       x <<= 2;
00275     }
00276   if (!(x & 0x80000000))
00277     {
00278       k++;
00279       if (!(x & 0x40000000))
00280         return 32;
00281     }
00282   return k;
00283 }
00284 
00285 int
00286 _DEFUN (lo0bits, (y), __ULong *y)
00287 {
00288   register int k;
00289   register __ULong x = *y;
00290 
00291   if (x & 7)
00292     {
00293       if (x & 1)
00294         return 0;
00295       if (x & 2)
00296         {
00297           *y = x >> 1;
00298           return 1;
00299         }
00300       *y = x >> 2;
00301       return 2;
00302     }
00303   k = 0;
00304   if (!(x & 0xffff))
00305     {
00306       k = 16;
00307       x >>= 16;
00308     }
00309   if (!(x & 0xff))
00310     {
00311       k += 8;
00312       x >>= 8;
00313     }
00314   if (!(x & 0xf))
00315     {
00316       k += 4;
00317       x >>= 4;
00318     }
00319   if (!(x & 0x3))
00320     {
00321       k += 2;
00322       x >>= 2;
00323     }
00324   if (!(x & 1))
00325     {
00326       k++;
00327       x >>= 1;
00328       if (!x & 1)
00329         return 32;
00330     }
00331   *y = x;
00332   return k;
00333 }
00334 
00335 _Bigint *
00336 _DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i)
00337 {
00338   _Bigint *b;
00339 
00340   b = Balloc (ptr, 1);
00341   b->_x[0] = i;
00342   b->_wds = 1;
00343   return b;
00344 }
00345 
00346 _Bigint *
00347 _DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b)
00348 {
00349   _Bigint *c;
00350   int k, wa, wb, wc;
00351   __ULong carry, y, z;
00352   __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00353 #ifdef Pack_32
00354   __ULong z2;
00355 #endif
00356 
00357   if (a->_wds < b->_wds)
00358     {
00359       c = a;
00360       a = b;
00361       b = c;
00362     }
00363   k = a->_k;
00364   wa = a->_wds;
00365   wb = b->_wds;
00366   wc = wa + wb;
00367   if (wc > a->_maxwds)
00368     k++;
00369   c = Balloc (ptr, k);
00370   for (x = c->_x, xa = x + wc; x < xa; x++)
00371     *x = 0;
00372   xa = a->_x;
00373   xae = xa + wa;
00374   xb = b->_x;
00375   xbe = xb + wb;
00376   xc0 = c->_x;
00377 #ifdef Pack_32
00378   for (; xb < xbe; xb++, xc0++)
00379     {
00380       if ((y = *xb & 0xffff) != 0)
00381         {
00382           x = xa;
00383           xc = xc0;
00384           carry = 0;
00385           do
00386             {
00387               z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
00388               carry = z >> 16;
00389               z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
00390               carry = z2 >> 16;
00391               Storeinc (xc, z2, z);
00392             }
00393           while (x < xae);
00394           *xc = carry;
00395         }
00396       if ((y = *xb >> 16) != 0)
00397         {
00398           x = xa;
00399           xc = xc0;
00400           carry = 0;
00401           z2 = *xc;
00402           do
00403             {
00404               z = (*x & 0xffff) * y + (*xc >> 16) + carry;
00405               carry = z >> 16;
00406               Storeinc (xc, z, z2);
00407               z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
00408               carry = z2 >> 16;
00409             }
00410           while (x < xae);
00411           *xc = z2;
00412         }
00413     }
00414 #else
00415   for (; xb < xbe; xc0++)
00416     {
00417       if (y = *xb++)
00418         {
00419           x = xa;
00420           xc = xc0;
00421           carry = 0;
00422           do
00423             {
00424               z = *x++ * y + *xc + carry;
00425               carry = z >> 16;
00426               *xc++ = z & 0xffff;
00427             }
00428           while (x < xae);
00429           *xc = carry;
00430         }
00431     }
00432 #endif
00433   for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
00434   c->_wds = wc;
00435   return c;
00436 }
00437 
00438 _Bigint *
00439 _DEFUN (pow5mult,
00440         (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
00441 {
00442   _Bigint *b1, *p5, *p51;
00443   int i;
00444   static _CONST int p05[3] = {5, 25, 125};
00445 
00446   if ((i = k & 3) != 0)
00447     b = multadd (ptr, b, p05[i - 1], 0);
00448 
00449   if (!(k >>= 2))
00450     return b;
00451   _REENT_CHECK_MP(ptr);
00452   if (!(p5 = _REENT_MP_P5S(ptr)))
00453     {
00454       /* first time */
00455       p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625);
00456       p5->_next = 0;
00457     }
00458   for (;;)
00459     {
00460       if (k & 1)
00461         {
00462           b1 = mult (ptr, b, p5);
00463           Bfree (ptr, b);
00464           b = b1;
00465         }
00466       if (!(k >>= 1))
00467         break;
00468       if (!(p51 = p5->_next))
00469         {
00470           p51 = p5->_next = mult (ptr, p5, p5);
00471           p51->_next = 0;
00472         }
00473       p5 = p51;
00474     }
00475   return b;
00476 }
00477 
00478 _Bigint *
00479 _DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k)
00480 {
00481   int i, k1, n, n1;
00482   _Bigint *b1;
00483   __ULong *x, *x1, *xe, z;
00484 
00485 #ifdef Pack_32
00486   n = k >> 5;
00487 #else
00488   n = k >> 4;
00489 #endif
00490   k1 = b->_k;
00491   n1 = n + b->_wds + 1;
00492   for (i = b->_maxwds; n1 > i; i <<= 1)
00493     k1++;
00494   b1 = Balloc (ptr, k1);
00495   x1 = b1->_x;
00496   for (i = 0; i < n; i++)
00497     *x1++ = 0;
00498   x = b->_x;
00499   xe = x + b->_wds;
00500 #ifdef Pack_32
00501   if (k &= 0x1f)
00502     {
00503       k1 = 32 - k;
00504       z = 0;
00505       do
00506         {
00507           *x1++ = *x << k | z;
00508           z = *x++ >> k1;
00509         }
00510       while (x < xe);
00511       if ((*x1 = z) != 0)
00512         ++n1;
00513     }
00514 #else
00515   if (k &= 0xf)
00516     {
00517       k1 = 16 - k;
00518       z = 0;
00519       do
00520         {
00521           *x1++ = *x << k & 0xffff | z;
00522           z = *x++ >> k1;
00523         }
00524       while (x < xe);
00525       if (*x1 = z)
00526         ++n1;
00527     }
00528 #endif
00529   else
00530     do
00531       *x1++ = *x++;
00532     while (x < xe);
00533   b1->_wds = n1 - 1;
00534   Bfree (ptr, b);
00535   return b1;
00536 }
00537 
00538 int
00539 _DEFUN (cmp, (a, b), _Bigint * a _AND _Bigint * b)
00540 {
00541   __ULong *xa, *xa0, *xb, *xb0;
00542   int i, j;
00543 
00544   i = a->_wds;
00545   j = b->_wds;
00546 #ifdef DEBUG
00547   if (i > 1 && !a->_x[i - 1])
00548     Bug ("cmp called with a->_x[a->_wds-1] == 0");
00549   if (j > 1 && !b->_x[j - 1])
00550     Bug ("cmp called with b->_x[b->_wds-1] == 0");
00551 #endif
00552   if (i -= j)
00553     return i;
00554   xa0 = a->_x;
00555   xa = xa0 + j;
00556   xb0 = b->_x;
00557   xb = xb0 + j;
00558   for (;;)
00559     {
00560       if (*--xa != *--xb)
00561         return *xa < *xb ? -1 : 1;
00562       if (xa <= xa0)
00563         break;
00564     }
00565   return 0;
00566 }
00567 
00568 _Bigint *
00569 _DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND
00570         _Bigint * a _AND _Bigint * b)
00571 {
00572   _Bigint *c;
00573   int i, wa, wb;
00574   __Long borrow, y;             /* We need signed shifts here. */
00575   __ULong *xa, *xae, *xb, *xbe, *xc;
00576 #ifdef Pack_32
00577   __Long z;
00578 #endif
00579 
00580   i = cmp (a, b);
00581   if (!i)
00582     {
00583       c = Balloc (ptr, 0);
00584       c->_wds = 1;
00585       c->_x[0] = 0;
00586       return c;
00587     }
00588   if (i < 0)
00589     {
00590       c = a;
00591       a = b;
00592       b = c;
00593       i = 1;
00594     }
00595   else
00596     i = 0;
00597   c = Balloc (ptr, a->_k);
00598   c->_sign = i;
00599   wa = a->_wds;
00600   xa = a->_x;
00601   xae = xa + wa;
00602   wb = b->_wds;
00603   xb = b->_x;
00604   xbe = xb + wb;
00605   xc = c->_x;
00606   borrow = 0;
00607 #ifdef Pack_32
00608   do
00609     {
00610       y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
00611       borrow = y >> 16;
00612       Sign_Extend (borrow, y);
00613       z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
00614       borrow = z >> 16;
00615       Sign_Extend (borrow, z);
00616       Storeinc (xc, z, y);
00617     }
00618   while (xb < xbe);
00619   while (xa < xae)
00620     {
00621       y = (*xa & 0xffff) + borrow;
00622       borrow = y >> 16;
00623       Sign_Extend (borrow, y);
00624       z = (*xa++ >> 16) + borrow;
00625       borrow = z >> 16;
00626       Sign_Extend (borrow, z);
00627       Storeinc (xc, z, y);
00628     }
00629 #else
00630   do
00631     {
00632       y = *xa++ - *xb++ + borrow;
00633       borrow = y >> 16;
00634       Sign_Extend (borrow, y);
00635       *xc++ = y & 0xffff;
00636     }
00637   while (xb < xbe);
00638   while (xa < xae)
00639     {
00640       y = *xa++ + borrow;
00641       borrow = y >> 16;
00642       Sign_Extend (borrow, y);
00643       *xc++ = y & 0xffff;
00644     }
00645 #endif
00646   while (!*--xc)
00647     wa--;
00648   c->_wds = wa;
00649   return c;
00650 }
00651 
00652 double
00653 _DEFUN (ulp, (_x), double _x)
00654 {
00655   union double_union x, a;
00656   register __Long L;
00657 
00658   x.d = _x;
00659 
00660   L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
00661 #ifndef Sudden_Underflow
00662   if (L > 0)
00663     {
00664 #endif
00665 #ifdef IBM
00666       L |= Exp_msk1 >> 4;
00667 #endif
00668       word0 (a) = L;
00669 #ifndef _DOUBLE_IS_32BITS
00670       word1 (a) = 0;
00671 #endif
00672 
00673 #ifndef Sudden_Underflow
00674     }
00675   else
00676     {
00677       L = -L >> Exp_shift;
00678       if (L < Exp_shift)
00679         {
00680           word0 (a) = 0x80000 >> L;
00681 #ifndef _DOUBLE_IS_32BITS
00682           word1 (a) = 0;
00683 #endif
00684         }
00685       else
00686         {
00687           word0 (a) = 0;
00688           L -= Exp_shift;
00689 #ifndef _DOUBLE_IS_32BITS
00690          word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
00691 #endif
00692         }
00693     }
00694 #endif
00695   return a.d;
00696 }
00697 
00698 double
00699 _DEFUN (b2d, (a, e),
00700         _Bigint * a _AND int *e)
00701 {
00702   __ULong *xa, *xa0, w, y, z;
00703   int k;
00704   union double_union d;
00705 #ifdef VAX
00706   __ULong d0, d1;
00707 #else
00708 #define d0 word0(d)
00709 #define d1 word1(d)
00710 #endif
00711 
00712   xa0 = a->_x;
00713   xa = xa0 + a->_wds;
00714   y = *--xa;
00715 #ifdef DEBUG
00716   if (!y)
00717     Bug ("zero y in b2d");
00718 #endif
00719   k = hi0bits (y);
00720   *e = 32 - k;
00721 #ifdef Pack_32
00722   if (k < Ebits)
00723     {
00724       d0 = Exp_1 | y >> (Ebits - k);
00725       w = xa > xa0 ? *--xa : 0;
00726 #ifndef _DOUBLE_IS_32BITS
00727       d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
00728 #endif
00729       goto ret_d;
00730     }
00731   z = xa > xa0 ? *--xa : 0;
00732   if (k -= Ebits)
00733     {
00734       d0 = Exp_1 | y << k | z >> (32 - k);
00735       y = xa > xa0 ? *--xa : 0;
00736 #ifndef _DOUBLE_IS_32BITS
00737       d1 = z << k | y >> (32 - k);
00738 #endif
00739     }
00740   else
00741     {
00742       d0 = Exp_1 | y;
00743 #ifndef _DOUBLE_IS_32BITS
00744       d1 = z;
00745 #endif
00746     }
00747 #else
00748   if (k < Ebits + 16)
00749     {
00750       z = xa > xa0 ? *--xa : 0;
00751       d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
00752       w = xa > xa0 ? *--xa : 0;
00753       y = xa > xa0 ? *--xa : 0;
00754       d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
00755       goto ret_d;
00756     }
00757   z = xa > xa0 ? *--xa : 0;
00758   w = xa > xa0 ? *--xa : 0;
00759   k -= Ebits + 16;
00760   d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
00761   y = xa > xa0 ? *--xa : 0;
00762   d1 = w << k + 16 | y << k;
00763 #endif
00764 ret_d:
00765 #ifdef VAX
00766   word0 (d) = d0 >> 16 | d0 << 16;
00767   word1 (d) = d1 >> 16 | d1 << 16;
00768 #else
00769 #undef d0
00770 #undef d1
00771 #endif
00772   return d.d;
00773 }
00774 
00775 _Bigint *
00776 _DEFUN (d2b,
00777         (ptr, _d, e, bits),
00778         struct _reent * ptr _AND
00779         double _d _AND
00780         int *e _AND
00781         int *bits)
00782 
00783 {
00784   union double_union d;
00785   _Bigint *b;
00786   int de, i, k;
00787   __ULong *x, y, z;
00788 #ifdef VAX
00789   __ULong d0, d1;
00790 #endif
00791   d.d = _d;
00792 #ifdef VAX
00793   d0 = word0 (d) >> 16 | word0 (d) << 16;
00794   d1 = word1 (d) >> 16 | word1 (d) << 16;
00795 #else
00796 #define d0 word0(d)
00797 #define d1 word1(d)
00798   d.d = _d;
00799 #endif
00800 
00801 #ifdef Pack_32
00802   b = Balloc (ptr, 1);
00803 #else
00804   b = Balloc (ptr, 2);
00805 #endif
00806   x = b->_x;
00807 
00808   z = d0 & Frac_mask;
00809   d0 &= 0x7fffffff;             /* clear sign bit, which we ignore */
00810 #ifdef Sudden_Underflow
00811   de = (int) (d0 >> Exp_shift);
00812 #ifndef IBM
00813   z |= Exp_msk11;
00814 #endif
00815 #else
00816   if ((de = (int) (d0 >> Exp_shift)) != 0)
00817     z |= Exp_msk1;
00818 #endif
00819 #ifdef Pack_32
00820 #ifndef _DOUBLE_IS_32BITS
00821   if (d1)
00822     {
00823       y = d1;
00824       k = lo0bits (&y);
00825       if (k)
00826         {
00827          x[0] = y | z << (32 - k);
00828           z >>= k;
00829         }
00830       else
00831         x[0] = y;
00832       i = b->_wds = (x[1] = z) ? 2 : 1;
00833     }
00834   else
00835 #endif
00836     {
00837 #ifdef DEBUG
00838       if (!z)
00839         Bug ("Zero passed to d2b");
00840 #endif
00841       k = lo0bits (&z);
00842       x[0] = z;
00843       i = b->_wds = 1;
00844 #ifndef _DOUBLE_IS_32BITS
00845       k += 32;
00846 #endif
00847     }
00848 #else
00849   if (d1)
00850     {
00851       y = d1;
00852       k = lo0bits (&y);
00853       if (k)
00854         if (k >= 16)
00855           {
00856             x[0] = y | z << 32 - k & 0xffff;
00857             x[1] = z >> k - 16 & 0xffff;
00858             x[2] = z >> k;
00859             i = 2;
00860           }
00861         else
00862           {
00863             x[0] = y & 0xffff;
00864             x[1] = y >> 16 | z << 16 - k & 0xffff;
00865             x[2] = z >> k & 0xffff;
00866             x[3] = z >> k + 16;
00867             i = 3;
00868           }
00869       else
00870         {
00871           x[0] = y & 0xffff;
00872           x[1] = y >> 16;
00873           x[2] = z & 0xffff;
00874           x[3] = z >> 16;
00875           i = 3;
00876         }
00877     }
00878   else
00879     {
00880 #ifdef DEBUG
00881       if (!z)
00882         Bug ("Zero passed to d2b");
00883 #endif
00884       k = lo0bits (&z);
00885       if (k >= 16)
00886         {
00887           x[0] = z;
00888           i = 0;
00889         }
00890       else
00891         {
00892           x[0] = z & 0xffff;
00893           x[1] = z >> 16;
00894           i = 1;
00895         }
00896       k += 32;
00897     }
00898   while (!x[i])
00899     --i;
00900   b->_wds = i + 1;
00901 #endif
00902 #ifndef Sudden_Underflow
00903   if (de)
00904     {
00905 #endif
00906 #ifdef IBM
00907       *e = (de - Bias - (P - 1) << 2) + k;
00908       *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
00909 #else
00910       *e = de - Bias - (P - 1) + k;
00911       *bits = P - k;
00912 #endif
00913 #ifndef Sudden_Underflow
00914     }
00915   else
00916     {
00917       *e = de - Bias - (P - 1) + 1 + k;
00918 #ifdef Pack_32
00919       *bits = 32 * i - hi0bits (x[i - 1]);
00920 #else
00921       *bits = (i + 2) * 16 - hi0bits (x[i]);
00922 #endif
00923     }
00924 #endif
00925   return b;
00926 }
00927 #undef d0
00928 #undef d1
00929 
00930 double
00931 _DEFUN (ratio, (a, b), _Bigint * a _AND _Bigint * b)
00932 
00933 {
00934   union double_union da, db;
00935   int k, ka, kb;
00936 
00937   da.d = b2d (a, &ka);
00938   db.d = b2d (b, &kb);
00939 #ifdef Pack_32
00940   k = ka - kb + 32 * (a->_wds - b->_wds);
00941 #else
00942   k = ka - kb + 16 * (a->_wds - b->_wds);
00943 #endif
00944 #ifdef IBM
00945   if (k > 0)
00946     {
00947       word0 (da) += (k >> 2) * Exp_msk1;
00948       if (k &= 3)
00949         da.d *= 1 << k;
00950     }
00951   else
00952     {
00953       k = -k;
00954       word0 (db) += (k >> 2) * Exp_msk1;
00955       if (k &= 3)
00956         db.d *= 1 << k;
00957     }
00958 #else
00959   if (k > 0)
00960     word0 (da) += k * Exp_msk1;
00961   else
00962     {
00963       k = -k;
00964       word0 (db) += k * Exp_msk1;
00965     }
00966 #endif
00967   return da.d / db.d;
00968 }
00969 
00970 
00971 _CONST double
00972   tens[] =
00973 {
00974   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
00975   1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
00976   1e20, 1e21, 1e22, 1e23, 1e24
00977 
00978 };
00979 
00980 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
00981 _CONST double bigtens[] =
00982 {1e16, 1e32, 1e64, 1e128, 1e256};
00983 
00984 _CONST double tinytens[] =
00985 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
00986 #else
00987 _CONST double bigtens[] =
00988 {1e16, 1e32};
00989 
00990 _CONST double tinytens[] =
00991 {1e-16, 1e-32};
00992 #endif
00993 
00994 
00995 double
00996 _DEFUN (_mprec_log10, (dig),
00997         int dig)
00998 {
00999   double v = 1.0;
01000   if (dig < 24)
01001     return tens[dig];
01002   while (dig > 0)
01003     {
01004       v *= 10;
01005       dig--;
01006     }
01007   return v;
01008 }
01009 
01010 #endif  // SMALL_LIB
01011 
01012 
01013 /*****************************************************************************************/
01014 /*                   FOR   SMALL_LIB                                                                                         */
01015 /*****************************************************************************************/
01016 
01017 
01018 #ifdef SMALL_LIB        
01019 
01020   /* 
01021         *       For SMALL_LIB, all references to balloc or bfree have been taken off 
01022    *    so as to avoided call of malloc libraries
01023         *       For each function returning a _Bigint * we added to the function a parameter
01024         *       of array type. The purpose of this array is to iniatialized pointers of _Bigint
01025         *       with the address of the array instead of the adress provided by balloc.
01026         *
01027         */
01028 
01029 
01030 
01031 _Bigint *
01032 _DEFUN (small_multadd, (ptr, b, m, a, tab),
01033         struct _reent *ptr _AND
01034         _Bigint * b _AND
01035         int m _AND
01036         int a _AND
01037         _Bigint tab[])
01038 {
01039    int i, wds;
01040   __ULong *x, y;
01041 #ifdef Pack_32
01042   __ULong xi, z;
01043 #endif
01044   _Bigint *b1;
01045 
01046   wds = b->_wds;
01047   x = b->_x;
01048   i = 0;
01049   do
01050     {
01051 #ifdef Pack_32
01052       xi = *x;
01053       y = (xi & 0xffff) * m + a;
01054       z = (xi >> 16) * m + (y >> 16);
01055       a = (int) (z >> 16);
01056       *x++ = (z << 16) + (y & 0xffff);
01057 #else
01058       y = *x * m + a;
01059       a = (int) (y >> 16);
01060       *x++ = y & 0xffff;
01061 #endif
01062     }
01063   while (++i < wds);
01064   if (a)
01065     {
01066       if (wds >= b->_maxwds)
01067         {
01068           b1=&tab[0];
01069           b1->_k=b->_k+1;
01070           b1->_maxwds = (1 <<(b->_k+1));
01071           b1->_sign=b1->_wds=0;
01072           Bcopy(b1,b);
01073           b=b1;
01074         }
01075       b->_x[wds++] = a;
01076       b->_wds = wds;
01077     }
01078   return b;
01079 }
01080 
01081 
01082 _Bigint *
01083 _DEFUN (small_s2b, (ptr, s, nd0, nd, y9,tab),
01084         struct _reent * ptr _AND
01085         _CONST char *s _AND
01086         int nd0 _AND
01087         int nd _AND
01088         __ULong y9 _AND
01089         _Bigint tab[])
01090 {
01091   _Bigint *b;
01092   _Bigint tab_b[50];
01093   int i, k;
01094   __Long x, y;
01095 
01096   x = (nd + 8) / 9;
01097   for (k = 0, y = 1; x > y; y <<= 1, k++);
01098 #ifdef Pack_32 
01099   b = &tab[0];
01100   b->_k=k;
01101   b->_maxwds = 1 <<k;
01102   b->_sign=0;
01103   b->_x[0] = y9;
01104   b->_wds = 1;
01105 #else
01106   b = &tab[0];
01107    b = &tab[0];
01108   b->_k=k+1;
01109   b->_mawxds = 1 <<(k+1);
01110   b->_sign=0;
01111   b->_x[0] = y9 & 0xffff;
01112   b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
01113 #endif
01114 
01115   i = 9;
01116   if (9 < nd0)
01117     {
01118       s += 9;
01119       
01120       while (++i < nd0);
01121       s++;
01122     }
01123   else
01124     s += 10;
01125   for (; i < nd; i++)
01126     b = small_multadd (ptr, b, 10, *s++ - '0',tab);
01127   return b;
01128 }
01129 
01130 int
01131 _DEFUN (small_hi0bits,
01132         (x), register __ULong x)
01133 {
01134   register int k = 0;
01135 
01136   if (!(x & 0xffff0000))
01137     {
01138       k = 16;
01139       x <<= 16;
01140     }
01141   if (!(x & 0xff000000))
01142     {
01143       k += 8;
01144       x <<= 8;
01145     }
01146   if (!(x & 0xf0000000))
01147     {
01148       k += 4;
01149       x <<= 4;
01150     }
01151   if (!(x & 0xc0000000))
01152     {
01153       k += 2;
01154       x <<= 2;
01155     }
01156   if (!(x & 0x80000000))
01157     {
01158       k++;
01159       if (!(x & 0x40000000))
01160         return 32;
01161     }
01162   return k;
01163 }
01164 
01165 int
01166 _DEFUN (small_lo0bits, (y), __ULong *y)
01167 {
01168   register int k;
01169   register __ULong x = *y;
01170 
01171   if (x & 7)
01172     {
01173       if (x & 1)
01174         return 0;
01175       if (x & 2)
01176         {
01177           *y = x >> 1;
01178           return 1;
01179         }
01180       *y = x >> 2;
01181       return 2;
01182     }
01183   k = 0;
01184   if (!(x & 0xffff))
01185     {
01186       k = 16;
01187       x >>= 16;
01188     }
01189   if (!(x & 0xff))
01190     {
01191       k += 8;
01192       x >>= 8;
01193     }
01194   if (!(x & 0xf))
01195     {
01196       k += 4;
01197       x >>= 4;
01198     }
01199   if (!(x & 0x3))
01200     {
01201       k += 2;
01202       x >>= 2;
01203     }
01204   if (!(x & 1))
01205     {
01206       k++;
01207       x >>= 1;
01208       if (!x & 1)
01209         return 32;
01210     }
01211   *y = x;
01212   return k;
01213 }
01214 
01215 
01216 _Bigint *
01217 _DEFUN (small_i2b, (ptr, i,tab), struct _reent * ptr _AND int i _AND _Bigint tab[])
01218 {
01219    _Bigint *b;
01220 
01221   b=&tab[0];
01222   b->_k=1;
01223   b->_maxwds = 1 << 1;
01224   b->_sign =0;
01225   b->_x[0] = i;
01226   b->_wds = 1;
01227   return b;
01228 }
01229 
01230 
01231 
01232 _Bigint *
01233 _DEFUN (small_mult, (ptr, a, b,tab), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b _AND _Bigint tab[])
01234 {
01235   _Bigint *c;
01236   int k, wa, wb, wc;
01237   __ULong carry, y, z;
01238   __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
01239 #ifdef Pack_32
01240   __ULong z2;
01241 #endif
01242 
01243   if (a->_wds < b->_wds)
01244     {
01245       c = a;
01246       a = b;
01247       b = c;
01248     }
01249   k = a->_k;
01250   wa = a->_wds;
01251   wb = b->_wds;
01252   wc = wa + wb;
01253   if (wc > a->_maxwds)
01254     k++;
01255   c=&tab[0];
01256   c->_k=k;
01257   c->_maxwds = 1 << k;
01258   c->_sign= c->_wds =0; 
01259   
01260   for (x = c->_x, xa = x + wc; x < xa; x++)
01261     *x = 0;
01262   xa = a->_x;
01263   xae = xa + wa;
01264   xb = b->_x;
01265   xbe = xb + wb;
01266   xc0 = c->_x;
01267 #ifdef Pack_32
01268   for (; xb < xbe; xb++, xc0++)
01269     {
01270       if ((y = *xb & 0xffff) != 0)
01271         {
01272           x = xa;
01273           xc = xc0;
01274           carry = 0;
01275           do
01276             {
01277               z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
01278               carry = z >> 16;
01279               z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
01280               carry = z2 >> 16;
01281               Storeinc (xc, z2, z);
01282             }
01283           while (x < xae);
01284           *xc = carry;
01285         }
01286       if ((y = *xb >> 16) != 0)
01287         {
01288           x = xa;
01289           xc = xc0;
01290           carry = 0;
01291           z2 = *xc;
01292           do
01293             {
01294               z = (*x & 0xffff) * y + (*xc >> 16) + carry;
01295               carry = z >> 16;
01296               Storeinc (xc, z, z2);
01297               z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
01298               carry = z2 >> 16;
01299             }
01300           while (x < xae);
01301           *xc = z2;
01302         }
01303     }
01304 #else
01305   for (; xb < xbe; xc0++)
01306     {
01307       if (y = *xb++)
01308         {
01309           x = xa;
01310           xc = xc0;
01311           carry = 0;
01312           do
01313             {
01314               z = *x++ * y + *xc + carry;
01315               carry = z >> 16;
01316               *xc++ = z & 0xffff;
01317             }
01318           while (x < xae);
01319           *xc = carry;
01320         }
01321     }
01322 #endif
01323   for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
01324   c->_wds = wc;
01325   return c;
01326 }
01327 
01328 
01329 
01330 _Bigint *
01331 _DEFUN (small_pow5mult,
01332         (ptr, b, k, tab), struct _reent * ptr _AND _Bigint * b _AND int k _AND _Bigint tab[])
01333 {
01334   _Bigint *b1, *p5, *p51;
01335   _Bigint tab_p5[50], tab_p5mult[50];
01336   _Bigint tab_b2[40], tab_b1[40];
01337   int i;
01338   int k_sauv;
01339   static _CONST int p05[3] = {5, 25, 125};
01340   if ((i = k & 3) != 0){
01341     k_sauv=k;
01342     k_sauv >>= 2;
01343     if (!(k_sauv)){ // we anticipate the case !(k>>=2) with return so as 
01344                                          //to provide the good address &tab[0]     
01345      b = small_multadd (ptr, b, p05[i - 1], 0,&tab[0]);
01346      return b;
01347     }
01348     else{
01349      b = small_multadd (ptr, b, p05[i - 1], 0,&tab_b1[0]); //&tab_b1[0] is an temporary address that
01350                                                                                                                                                   // that we provide to b
01351     }
01352   }   
01353   if (!(k>>= 2)){
01354     return b; 
01355   }   
01356   if (p5 != tab_p5 ){ /* first time */
01357       p5 =small_i2b (ptr, 625,tab_p5); 
01358       p5->_next = 0;
01359   }
01360    
01361   for (;;) {// We are in a loop for each variable and call of mult we must provide a
01362                            // an address which differs from the current address of the variable
01363                            // that is why we add a test on the current address
01364 
01365   if (k & 1){     
01366           
01367           k_sauv = k ;
01368           if (!(k_sauv >>= 1) ){ // it is the last passage in the loop
01369                                                                          // we must provide the good address
01370             b1 = small_mult (ptr, b, p5,&tab[0]);
01371             b  = b1;      
01372             break;
01373           }
01374           else {
01375             if (b == &tab_b1[0]){
01376                b1 = small_mult (ptr, b, p5,&tab_b2[0]);
01377             }
01378             else {
01379                b1 = small_mult (ptr, b, p5,&tab_b1[0]);
01380             }  
01381           }
01382           b  = b1;                
01383   }
01384   if (!(k >>= 1))
01385         break;
01386       if (!(p51 = p5->_next))
01387         {
01388           if ( p5 == &tab_p5[0] ){
01389           p51 = p5->_next = small_mult (ptr, p5, p5,&tab_p5mult[0]);
01390           }
01391           else {  
01392           p51 = p5->_next = small_mult (ptr, p5, p5,&tab_p5[0]);
01393           }       
01394           p51->_next = 0;
01395         }
01396       p5 = p51;
01397     }
01398   return b;
01399 }
01400 
01401 
01402 _Bigint *
01403 _DEFUN (small_lshift, (ptr, b, k,tab), struct _reent * ptr _AND _Bigint * b _AND int k _AND _Bigint tab[] )
01404 {
01405   int i, k1, n, n1;
01406   _Bigint *b1;
01407   __ULong *x, *x1, *xe, z;
01408 
01409 #ifdef Pack_32
01410   n = k >> 5;
01411 #else
01412   n = k >> 4;
01413 #endif
01414   k1 = b->_k;
01415   n1 = n + b->_wds + 1;
01416   for (i = b->_maxwds; n1 > i; i <<= 1)
01417   k1++;
01418  //b1 = Balloc (ptr, k1);
01419  b1=&tab[0];
01420  b1->_k=k1;
01421  b1->_maxwds = 1 << k1;
01422  b1->_sign= b1->_wds =0;
01423   x1 = b1->_x;
01424   for (i = 0; i < n; i++)
01425     *x1++ = 0;
01426   x = b->_x;
01427   xe = x + b->_wds;
01428 #ifdef Pack_32
01429   if (k &= 0x1f)
01430     {
01431       k1 = 32 - k;
01432       z = 0;
01433       do
01434         {
01435           *x1++ = *x << k | z;
01436           z = *x++ >> k1;
01437         }
01438       while (x < xe);
01439       if ((*x1 = z) != 0)
01440         ++n1;
01441     }
01442 #else
01443   if (k &= 0xf)
01444     {
01445       k1 = 16 - k;
01446       z = 0;
01447       do
01448         {
01449           *x1++ = *x << k & 0xffff | z;
01450           z = *x++ >> k1;
01451         }
01452       while (x < xe);
01453       if (*x1 = z)
01454         ++n1;
01455     }
01456 #endif
01457   else
01458     do
01459       *x1++ = *x++;
01460     while (x < xe);
01461   b1->_wds = n1 - 1;
01462   //Bfree (ptr, b);
01463   return b1;
01464 }
01465 
01466 int
01467 _DEFUN (small_cmp, (a, b), _Bigint * a _AND _Bigint * b)
01468 {
01469   __ULong *xa, *xa0, *xb, *xb0;
01470   int i, j;
01471 
01472   i = a->_wds;
01473   j = b->_wds;
01474 #ifdef DEBUG
01475   if (i > 1 && !a->_x[i - 1])
01476     Bug ("cmp called with a->_x[a->_wds-1] == 0");
01477   if (j > 1 && !b->_x[j - 1])
01478     Bug ("cmp called with b->_x[b->_wds-1] == 0");
01479 #endif
01480   if (i -= j)
01481     return i;
01482   xa0 = a->_x;
01483   xa = xa0 + j;
01484   xb0 = b->_x;
01485   xb = xb0 + j;
01486   for (;;)
01487     {
01488       if (*--xa != *--xb)
01489         return *xa < *xb ? -1 : 1;
01490       if (xa <= xa0)
01491         break;
01492     }
01493   return 0;
01494 }
01495 
01496 
01497 
01498 _Bigint *
01499 _DEFUN (small_diff, (ptr, a, b,tab), struct _reent * ptr _AND
01500         _Bigint * a _AND _Bigint * b _AND _Bigint tab[])
01501 {
01502   _Bigint *c;
01503   int i, wa, wb;
01504   __Long borrow, y;             /* We need signed shifts here. */
01505   __ULong *xa, *xae, *xb, *xbe, *xc;
01506 #ifdef Pack_32
01507   __Long z;
01508 #endif
01509 
01510   i = small_cmp (a, b);
01511   if (!i)
01512     {
01513       c=&tab[0];
01514       c->_k = 0;
01515       c->_maxwds = (1 << 0) ;
01516       c->_sign = 0 ;
01517       c->_wds = 1;
01518       c->_x[0] = 0;
01519       return c;
01520     }
01521   if (i < 0)
01522     {
01523       c = a;
01524       a = b;
01525       b = c;
01526       i = 1;
01527     }
01528   else
01529   i = 0;
01530   c=&tab[0];
01531   c->_k = a->_k;
01532   c->_maxwds = (1 << (a->_k)) ;
01533   c->_wds = 0 ;
01534   c->_sign = i;
01535   wa = a->_wds;
01536   xa = a->_x;
01537   xae = xa + wa;
01538   wb = b->_wds;
01539   xb = b->_x;
01540   xbe = xb + wb;
01541   xc = c->_x;
01542   borrow = 0;
01543 #ifdef Pack_32
01544   do
01545     {
01546       y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
01547       borrow = y >> 16;
01548       Sign_Extend (borrow, y);
01549       z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
01550       borrow = z >> 16;
01551       Sign_Extend (borrow, z);
01552       Storeinc (xc, z, y);
01553     }
01554   while (xb < xbe);
01555   while (xa < xae)
01556     {
01557       y = (*xa & 0xffff) + borrow;
01558       borrow = y >> 16;
01559       Sign_Extend (borrow, y);
01560       z = (*xa++ >> 16) + borrow;
01561       borrow = z >> 16;
01562       Sign_Extend (borrow, z);
01563       Storeinc (xc, z, y);
01564     }
01565 #else
01566   do
01567     {
01568       y = *xa++ - *xb++ + borrow;
01569       borrow = y >> 16;
01570       Sign_Extend (borrow, y);
01571       *xc++ = y & 0xffff;
01572     }
01573   while (xb < xbe);
01574   while (xa < xae)
01575     {
01576       y = *xa++ + borrow;
01577       borrow = y >> 16;
01578       Sign_Extend (borrow, y);
01579       *xc++ = y & 0xffff;
01580     }
01581 #endif
01582   while (!*--xc)
01583     wa--;
01584   c->_wds = wa;
01585   return c;
01586 }
01587 
01588 double
01589 _DEFUN (small_ulp, (_x), double _x)
01590 {
01591   union double_union x, a;
01592   register __Long L;
01593 
01594   x.d = _x;
01595 
01596   L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
01597 #ifndef Sudden_Underflow
01598   if (L > 0)
01599     {
01600 #endif
01601 #ifdef IBM
01602       L |= Exp_msk1 >> 4;
01603 #endif
01604       word0 (a) = L;
01605 #ifndef _DOUBLE_IS_32BITS
01606       word1 (a) = 0;
01607 #endif
01608 
01609 #ifndef Sudden_Underflow
01610     }
01611   else
01612     {
01613       L = -L >> Exp_shift;
01614       if (L < Exp_shift)
01615         {
01616           word0 (a) = 0x80000 >> L;
01617 #ifndef _DOUBLE_IS_32BITS
01618           word1 (a) = 0;
01619 #endif
01620         }
01621       else
01622         {
01623           word0 (a) = 0;
01624           L -= Exp_shift;
01625 #ifndef _DOUBLE_IS_32BITS
01626          word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
01627 #endif
01628         }
01629     }
01630 #endif
01631   return a.d;
01632 }
01633 
01634 double
01635 _DEFUN (small_b2d, (a, e),
01636         _Bigint * a _AND int *e)
01637 {
01638   __ULong *xa, *xa0, w, y, z;
01639   int k;
01640   union double_union d;
01641 #ifdef VAX
01642   __ULong d0, d1;
01643 #else
01644 #define d0 word0(d)
01645 #define d1 word1(d)
01646 #endif
01647 
01648   xa0 = a->_x;
01649   xa = xa0 + a->_wds;
01650   y = *--xa;
01651 #ifdef DEBUG
01652   if (!y)
01653     Bug ("zero y in b2d");
01654 #endif
01655   k = small_hi0bits (y);
01656   *e = 32 - k;
01657 #ifdef Pack_32
01658   if (k < Ebits)
01659     {
01660       d0 = Exp_1 | y >> (Ebits - k);
01661       w = xa > xa0 ? *--xa : 0;
01662 #ifndef _DOUBLE_IS_32BITS
01663       d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
01664 #endif
01665       goto ret_d;
01666     }
01667   z = xa > xa0 ? *--xa : 0;
01668   if (k -= Ebits)
01669     {
01670       d0 = Exp_1 | y << k | z >> (32 - k);
01671       y = xa > xa0 ? *--xa : 0;
01672 #ifndef _DOUBLE_IS_32BITS
01673       d1 = z << k | y >> (32 - k);
01674 #endif
01675     }
01676   else
01677     {
01678       d0 = Exp_1 | y;
01679 #ifndef _DOUBLE_IS_32BITS
01680       d1 = z;
01681 #endif
01682     }
01683 #else
01684   if (k < Ebits + 16)
01685     {
01686       z = xa > xa0 ? *--xa : 0;
01687       d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
01688       w = xa > xa0 ? *--xa : 0;
01689       y = xa > xa0 ? *--xa : 0;
01690       d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
01691       goto ret_d;
01692     }
01693   z = xa > xa0 ? *--xa : 0;
01694   w = xa > xa0 ? *--xa : 0;
01695   k -= Ebits + 16;
01696   d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
01697   y = xa > xa0 ? *--xa : 0;
01698   d1 = w << k + 16 | y << k;
01699 #endif
01700 ret_d:
01701 #ifdef VAX
01702   word0 (d) = d0 >> 16 | d0 << 16;
01703   word1 (d) = d1 >> 16 | d1 << 16;
01704 #else
01705 #undef d0
01706 #undef d1
01707 #endif
01708   return d.d;
01709 }
01710 
01711 
01712 _Bigint *
01713 _DEFUN (small_d2b,
01714         (ptr, _d, e, bits,tab),
01715         struct _reent * ptr _AND
01716         double _d _AND
01717         int *e _AND
01718         int *bits _AND
01719         _Bigint tab[]
01720         )
01721 
01722 {
01723   union double_union d;
01724   _Bigint *b;
01725   
01726   int de, i, k;
01727   __ULong *x, y, z;
01728 #ifdef VAX
01729   __ULong d0, d1;
01730 #endif
01731   d.d = _d;
01732 #ifdef VAX
01733   d0 = word0 (d) >> 16 | word0 (d) << 16;
01734   d1 = word1 (d) >> 16 | word1 (d) << 16;
01735 #else
01736 #define d0 word0(d)
01737 #define d1 word1(d)
01738   d.d = _d;
01739 #endif
01740 
01741 #ifdef Pack_32
01742   b=&tab[0];
01743   b->_k = 1;
01744   b->_maxwds = (1 << 1) ;
01745   b->_sign = b->_wds = 0 ;
01746   
01747 #else
01748   b=&tab[0];
01749   b->_k = 2;
01750   b->_maxwds = (1 << 2) ;
01751   b->_sign = b->_wds = 0 ;
01752 #endif 
01753 
01754 
01755   x = b->_x;
01756 
01757   z = d0 & Frac_mask;
01758   d0 &= 0x7fffffff;             /* clear sign bit, which we ignore */
01759 #ifdef Sudden_Underflow
01760   de = (int) (d0 >> Exp_shift);
01761 #ifndef IBM
01762   z |= Exp_msk11;
01763 #endif
01764 #else
01765   if ((de = (int) (d0 >> Exp_shift)) != 0)
01766     z |= Exp_msk1;
01767 #endif
01768 #ifdef Pack_32
01769 #ifndef _DOUBLE_IS_32BITS
01770   if (d1)
01771     {
01772       y = d1;
01773       k = small_lo0bits (&y);
01774       if (k)
01775         {
01776          x[0] = y | z << (32 - k);
01777           z >>= k;
01778         }
01779       else
01780         x[0] = y;
01781       i = b->_wds = (x[1] = z) ? 2 : 1;
01782     }
01783   else
01784 #endif
01785     {
01786 #ifdef DEBUG
01787       if (!z)
01788         Bug ("Zero passed to d2b");
01789 #endif
01790       k = small_lo0bits (&z);
01791       x[0] = z;
01792       i = b->_wds = 1;
01793 #ifndef _DOUBLE_IS_32BITS
01794       k += 32;
01795 #endif
01796     }
01797 #else
01798   if (d1)
01799     {
01800       y = d1;
01801       k = small_lo0bits (&y);
01802       if (k)
01803         if (k >= 16)
01804           {
01805             x[0] = y | z << 32 - k & 0xffff;
01806             x[1] = z >> k - 16 & 0xffff;
01807             x[2] = z >> k;
01808             i = 2;
01809           }
01810         else
01811           {
01812             x[0] = y & 0xffff;
01813             x[1] = y >> 16 | z << 16 - k & 0xffff;
01814             x[2] = z >> k & 0xffff;
01815             x[3] = z >> k + 16;
01816             i = 3;
01817           }
01818       else
01819         {
01820           x[0] = y & 0xffff;
01821           x[1] = y >> 16;
01822           x[2] = z & 0xffff;
01823           x[3] = z >> 16;
01824           i = 3;
01825         }
01826     }
01827   else
01828     {
01829 #ifdef DEBUG
01830       if (!z)
01831         Bug ("Zero passed to d2b");
01832 #endif
01833       k = lo0bits (&z);
01834       if (k >= 16)
01835         {
01836           x[0] = z;
01837           i = 0;
01838         }
01839       else
01840         {
01841           x[0] = z & 0xffff;
01842           x[1] = z >> 16;
01843           i = 1;
01844         }
01845       k += 32;
01846     }
01847   while (!x[i])
01848     --i;
01849   b->_wds = i + 1;
01850 #endif
01851 #ifndef Sudden_Underflow
01852   if (de)
01853     {
01854 #endif
01855 #ifdef IBM
01856       *e = (de - Bias - (P - 1) << 2) + k;
01857       *bits = 4 * P + 8 - k - small_hi0bits (word0 (d) & Frac_mask);
01858 #else
01859       *e = de - Bias - (P - 1) + k;
01860       *bits = P - k;
01861 #endif
01862 #ifndef Sudden_Underflow
01863     }
01864   else
01865     {
01866       *e = de - Bias - (P - 1) + 1 + k;
01867 #ifdef Pack_32
01868       *bits = 32 * i - small_hi0bits (x[i - 1]);
01869 #else
01870       *bits = (i + 2) * 16 - small_hi0bits (x[i]);
01871 #endif
01872     }
01873 #endif
01874   return b;
01875 }
01876 #undef d0
01877 #undef d1
01878 
01879 double
01880 _DEFUN (small_ratio, (a, b), _Bigint * a _AND _Bigint * b)
01881 
01882 {
01883   union double_union da, db;
01884   int k, ka, kb;
01885 
01886   da.d = small_b2d (a, &ka);
01887   db.d = small_b2d (b, &kb);
01888 #ifdef Pack_32
01889   k = ka - kb + 32 * (a->_wds - b->_wds);
01890 #else
01891   k = ka - kb + 16 * (a->_wds - b->_wds);
01892 #endif
01893 #ifdef IBM
01894   if (k > 0)
01895     {
01896       word0 (da) += (k >> 2) * Exp_msk1;
01897       if (k &= 3)
01898         da.d *= 1 << k;
01899     }
01900   else
01901     {
01902       k = -k;
01903       word0 (db) += (k >> 2) * Exp_msk1;
01904       if (k &= 3)
01905         db.d *= 1 << k;
01906     }
01907 #else
01908   if (k > 0)
01909     word0 (da) += k * Exp_msk1;
01910   else
01911     {
01912       k = -k;
01913       word0 (db) += k * Exp_msk1;
01914     }
01915 #endif
01916   return da.d / db.d;
01917 }
01918 
01919 
01920 _CONST double
01921   small_tens[] =
01922 {
01923   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01924   1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01925   1e20, 1e21, 1e22, 1e23, 1e24
01926 
01927 };
01928 
01929 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
01930 _CONST double small_bigtens[] =
01931 {1e16, 1e32, 1e64, 1e128, 1e256};
01932 
01933 _CONST double small_tinytens[] =
01934 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
01935 #else
01936 _CONST double small_bigtens[] =
01937 {1e16, 1e32};
01938 
01939 _CONST double small_tinytens[] =
01940 {1e-16, 1e-32};
01941 #endif
01942 
01943 
01944 
01945 
01946 double
01947 _DEFUN (small__mprec_log10, (dig),
01948         int dig)
01949 {
01950   double v = 1.0;
01951   if (dig < 24)
01952     return small_tens[dig];
01953   while (dig > 0)
01954     {
01955       v *= 10;
01956       dig--;
01957     }
01958   return v;
01959 }
01960 
01961 #endif // SMALL_PRINTF
01962 
01963