newlib/libc/stdlib/mprec.c

   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 #ifdef _REENT_THREAD_LOCAL
  90 _Thread_local struct _Bigint *_tls_mp_p5s;
  91 _Thread_local struct _Bigint **_tls_mp_freelist;
  92 #endif
  93
  94 /* This is defined in sys/reent.h as (sizeof (size_t) << 3) now, as in NetBSD.
  95    The old value of 15 was wrong and made newlib vulnerable against buffer
  96    overrun attacks (CVE-2009-0689), same as other implementations of gdtoa
  97    based on BSD code.
  98 #define _Kmax 15
  99 */
 100
 101 _Bigint *
 102 Balloc (struct _reent *ptr, int k)
 103 {
 104   int x;
 105   _Bigint *rv ;
 106
 107   _REENT_CHECK_MP(ptr);
 108   if (_REENT_MP_FREELIST(ptr) == NULL)
 109     {
 110       /* Allocate a list of pointers to the mprec objects */
 111       _REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr,
 112                                                       sizeof (struct _Bigint *),
 113                                                       _Kmax + 1);
 114       if (_REENT_MP_FREELIST(ptr) == NULL)
 115         {
 116           return NULL;
 117         }
 118     }
 119
 120   if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0)
 121     {
 122       _REENT_MP_FREELIST(ptr)[k] = rv->_next;
 123     }
 124   else
 125     {
 126       x = 1 << k;
 127       /* Allocate an mprec Bigint and stick in in the freelist */
 128       rv = (_Bigint *) _calloc_r (ptr,
 129                                   1,
 130                                   sizeof (_Bigint) +
 131                                   (x-1) * sizeof(rv->_x));
 132       if (rv == NULL) return NULL;
 133       rv->_k = k;
 134       rv->_maxwds = x;
 135     }
 136   rv->_sign = rv->_wds = 0;
 137   return rv;
 138 }
 139
 140 void
 141 Bfree (struct _reent *ptr, _Bigint * v)
 142 {
 143   _REENT_CHECK_MP(ptr);
 144   if (v)
 145     {
 146       v->_next = _REENT_MP_FREELIST(ptr)[v->_k];
 147       _REENT_MP_FREELIST(ptr)[v->_k] = v;
 148     }
 149 }
 150
 151 _Bigint *
 152 multadd (struct _reent *ptr,
 153         _Bigint * b,
 154         int m,
 155         int a)
 156 {
 157   int i, wds;
 158   __ULong *x, y;
 159 #ifdef Pack_32
 160   __ULong xi, z;
 161 #endif
 162   _Bigint *b1;
 163
 164   wds = b->_wds;
 165   x = b->_x;
 166   i = 0;
 167   do
 168     {
 169 #ifdef Pack_32
 170       xi = *x;
 171       y = (xi & 0xffff) * m + a;
 172       z = (xi >> 16) * m + (y >> 16);
 173       a = (int) (z >> 16);
 174       *x++ = (z << 16) + (y & 0xffff);
 175 #else
 176       y = *x * m + a;
 177       a = (int) (y >> 16);
 178       *x++ = y & 0xffff;
 179 #endif
 180     }
 181   while (++i < wds);
 182   if (a)
 183     {
 184       if (wds >= b->_maxwds)
 185         {
 186           b1 = eBalloc (ptr, b->_k + 1);
 187           Bcopy (b1, b);
 188           Bfree (ptr, b);
 189           b = b1;
 190         }
 191       b->_x[wds++] = a;
 192       b->_wds = wds;
 193     }
 194   return b;
 195 }
 196
 197 _Bigint *
 198 s2b (struct _reent * ptr,
 199         const char *s,
 200         int nd0,
 201         int nd,
 202         __ULong y9)
 203 {
 204   _Bigint *b;
 205   int i, k;
 206   __Long x, y;
 207
 208   x = (nd + 8) / 9;
 209   for (k = 0, y = 1; x > y; y <<= 1, k++);
 210 #ifdef Pack_32
 211   b = eBalloc (ptr, k);
 212   b->_x[0] = y9;
 213   b->_wds = 1;
 214 #else
 215   b = eBalloc (ptr, k + 1);
 216   b->_x[0] = y9 & 0xffff;
 217   b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
 218 #endif
 219
 220   i = 9;
 221   if (9 < nd0)
 222     {
 223       s += 9;
 224       do
 225         b = multadd (ptr, b, 10, *s++ - '0');
 226       while (++i < nd0);
 227       s++;
 228     }
 229   else
 230     s += 10;
 231   for (; i < nd; i++)
 232     b = multadd (ptr, b, 10, *s++ - '0');
 233   return b;
 234 }
 235
 236 int
 237 hi0bits (register __ULong x)
 238 {
 239   register int k = 0;
 240
 241   if (!(x & 0xffff0000))
 242     {
 243       k = 16;
 244       x <<= 16;
 245     }
 246   if (!(x & 0xff000000))
 247     {
 248       k += 8;
 249       x <<= 8;
 250     }
 251   if (!(x & 0xf0000000))
 252     {
 253       k += 4;
 254       x <<= 4;
 255     }
 256   if (!(x & 0xc0000000))
 257     {
 258       k += 2;
 259       x <<= 2;
 260     }
 261   if (!(x & 0x80000000))
 262     {
 263       k++;
 264       if (!(x & 0x40000000))
 265         return 32;
 266     }
 267   return k;
 268 }
 269
 270 int
 271 lo0bits (__ULong *y)
 272 {
 273   register int k;
 274   register __ULong x = *y;
 275
 276   if (x & 7)
 277     {
 278       if (x & 1)
 279         return 0;
 280       if (x & 2)
 281         {
 282           *y = x >> 1;
 283           return 1;
 284         }
 285       *y = x >> 2;
 286       return 2;
 287     }
 288   k = 0;
 289   if (!(x & 0xffff))
 290     {
 291       k = 16;
 292       x >>= 16;
 293     }
 294   if (!(x & 0xff))
 295     {
 296       k += 8;
 297       x >>= 8;
 298     }
 299   if (!(x & 0xf))
 300     {
 301       k += 4;
 302       x >>= 4;
 303     }
 304   if (!(x & 0x3))
 305     {
 306       k += 2;
 307       x >>= 2;
 308     }
 309   if (!(x & 1))
 310     {
 311       k++;
 312       x >>= 1;
 313       if (!x & 1)
 314         return 32;
 315     }
 316   *y = x;
 317   return k;
 318 }
 319
 320 _Bigint *
 321 i2b (struct _reent * ptr, int i)
 322 {
 323   _Bigint *b;
 324
 325   b = eBalloc (ptr, 1);
 326   b->_x[0] = i;
 327   b->_wds = 1;
 328   return b;
 329 }
 330
 331 _Bigint *
 332 mult (struct _reent * ptr, _Bigint * a, _Bigint * b)
 333 {
 334   _Bigint *c;
 335   int k, wa, wb, wc;
 336   __ULong carry, y, z;
 337   __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 338 #ifdef Pack_32
 339   __ULong z2;
 340 #endif
 341
 342   if (a->_wds < b->_wds)
 343     {
 344       c = a;
 345       a = b;
 346       b = c;
 347     }
 348   k = a->_k;
 349   wa = a->_wds;
 350   wb = b->_wds;
 351   wc = wa + wb;
 352   if (wc > a->_maxwds)
 353     k++;
 354   c = eBalloc (ptr, k);
 355   for (x = c->_x, xa = x + wc; x < xa; x++)
 356     *x = 0;
 357   xa = a->_x;
 358   xae = xa + wa;
 359   xb = b->_x;
 360   xbe = xb + wb;
 361   xc0 = c->_x;
 362 #ifdef Pack_32
 363   for (; xb < xbe; xb++, xc0++)
 364     {
 365       if ((y = *xb & 0xffff) != 0)
 366         {
 367           x = xa;
 368           xc = xc0;
 369           carry = 0;
 370           do
 371             {
 372               z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 373               carry = z >> 16;
 374               z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 375               carry = z2 >> 16;
 376               Storeinc (xc, z2, z);
 377             }
 378           while (x < xae);
 379           *xc = carry;
 380         }
 381       if ((y = *xb >> 16) != 0)
 382         {
 383           x = xa;
 384           xc = xc0;
 385           carry = 0;
 386           z2 = *xc;
 387           do
 388             {
 389               z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 390               carry = z >> 16;
 391               Storeinc (xc, z, z2);
 392               z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 393               carry = z2 >> 16;
 394             }
 395           while (x < xae);
 396           *xc = z2;
 397         }
 398     }
 399 #else
 400   for (; xb < xbe; xc0++)
 401     {
 402       if (y = *xb++)
 403         {
 404           x = xa;
 405           xc = xc0;
 406           carry = 0;
 407           do
 408             {
 409               z = *x++ * y + *xc + carry;
 410               carry = z >> 16;
 411               *xc++ = z & 0xffff;
 412             }
 413           while (x < xae);
 414           *xc = carry;
 415         }
 416     }
 417 #endif
 418   for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
 419   c->_wds = wc;
 420   return c;
 421 }
 422
 423 _Bigint *
 424 pow5mult (struct _reent * ptr, _Bigint * b, int k)
 425 {
 426   _Bigint *b1, *p5, *p51;
 427   int i;
 428   static const int p05[3] = {5, 25, 125};
 429
 430   if ((i = k & 3) != 0)
 431     b = multadd (ptr, b, p05[i - 1], 0);
 432
 433   if (!(k >>= 2))
 434     return b;
 435   _REENT_CHECK_MP(ptr);
 436   if (!(p5 = _REENT_MP_P5S(ptr)))
 437     {
 438       /* first time */
 439       p5 = _REENT_MP_P5S(ptr) = i2b (ptr, 625);
 440       p5->_next = 0;
 441     }
 442   for (;;)
 443     {
 444       if (k & 1)
 445         {
 446           b1 = mult (ptr, b, p5);
 447           Bfree (ptr, b);
 448           b = b1;
 449         }
 450       if (!(k >>= 1))
 451         break;
 452       if (!(p51 = p5->_next))
 453         {
 454           p51 = p5->_next = mult (ptr, p5, p5);
 455           p51->_next = 0;
 456         }
 457       p5 = p51;
 458     }
 459   return b;
 460 }
 461
 462 _Bigint *
 463 lshift (struct _reent * ptr, _Bigint * b, int k)
 464 {
 465   int i, k1, n, n1;
 466   _Bigint *b1;
 467   __ULong *x, *x1, *xe, z;
 468
 469 #ifdef Pack_32
 470   n = k >> 5;
 471 #else
 472   n = k >> 4;
 473 #endif
 474   k1 = b->_k;
 475   n1 = n + b->_wds + 1;
 476   for (i = b->_maxwds; n1 > i; i <<= 1)
 477     k1++;
 478   b1 = eBalloc (ptr, k1);
 479   x1 = b1->_x;
 480   for (i = 0; i < n; i++)
 481     *x1++ = 0;
 482   x = b->_x;
 483   xe = x + b->_wds;
 484 #ifdef Pack_32
 485   if (k &= 0x1f)
 486     {
 487       k1 = 32 - k;
 488       z = 0;
 489       do
 490         {
 491           *x1++ = *x << k | z;
 492           z = *x++ >> k1;
 493         }
 494       while (x < xe);
 495       if ((*x1 = z) != 0)
 496         ++n1;
 497     }
 498 #else
 499   if (k &= 0xf)
 500     {
 501       k1 = 16 - k;
 502       z = 0;
 503       do
 504         {
 505           *x1++ = *x << k & 0xffff | z;
 506           z = *x++ >> k1;
 507         }
 508       while (x < xe);
 509       if (*x1 = z)
 510         ++n1;
 511     }
 512 #endif
 513   else
 514     do
 515       *x1++ = *x++;
 516     while (x < xe);
 517   b1->_wds = n1 - 1;
 518   Bfree (ptr, b);
 519   return b1;
 520 }
 521
 522 int
 523 cmp (_Bigint * a, _Bigint * b)
 524 {
 525   __ULong *xa, *xa0, *xb, *xb0;
 526   int i, j;
 527
 528   i = a->_wds;
 529   j = b->_wds;
 530 #ifdef DEBUG
 531   if (i > 1 && !a->_x[i - 1])
 532     Bug ("cmp called with a->_x[a->_wds-1] == 0");
 533   if (j > 1 && !b->_x[j - 1])
 534     Bug ("cmp called with b->_x[b->_wds-1] == 0");
 535 #endif
 536   if (i -= j)
 537     return i;
 538   xa0 = a->_x;
 539   xa = xa0 + j;
 540   xb0 = b->_x;
 541   xb = xb0 + j;
 542   for (;;)
 543     {
 544       if (*--xa != *--xb)
 545         return *xa < *xb ? -1 : 1;
 546       if (xa <= xa0)
 547         break;
 548     }
 549   return 0;
 550 }
 551
 552 _Bigint *
 553 diff (struct _reent * ptr,
 554         _Bigint * a, _Bigint * b)
 555 {
 556   _Bigint *c;
 557   int i, wa, wb;
 558   __Long borrow, y;             /* We need signed shifts here. */
 559   __ULong *xa, *xae, *xb, *xbe, *xc;
 560 #ifdef Pack_32
 561   __Long z;
 562 #endif
 563
 564   i = cmp (a, b);
 565   if (!i)
 566     {
 567       c = eBalloc (ptr, 0);
 568       c->_wds = 1;
 569       c->_x[0] = 0;
 570       return c;
 571     }
 572   if (i < 0)
 573     {
 574       c = a;
 575       a = b;
 576       b = c;
 577       i = 1;
 578     }
 579   else
 580     i = 0;
 581   c = eBalloc (ptr, a->_k);
 582   c->_sign = i;
 583   wa = a->_wds;
 584   xa = a->_x;
 585   xae = xa + wa;
 586   wb = b->_wds;
 587   xb = b->_x;
 588   xbe = xb + wb;
 589   xc = c->_x;
 590   borrow = 0;
 591 #ifdef Pack_32
 592   do
 593     {
 594       y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
 595       borrow = y >> 16;
 596       Sign_Extend (borrow, y);
 597       z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
 598       borrow = z >> 16;
 599       Sign_Extend (borrow, z);
 600       Storeinc (xc, z, y);
 601     }
 602   while (xb < xbe);
 603   while (xa < xae)
 604     {
 605       y = (*xa & 0xffff) + borrow;
 606       borrow = y >> 16;
 607       Sign_Extend (borrow, y);
 608       z = (*xa++ >> 16) + borrow;
 609       borrow = z >> 16;
 610       Sign_Extend (borrow, z);
 611       Storeinc (xc, z, y);
 612     }
 613 #else
 614   do
 615     {
 616       y = *xa++ - *xb++ + borrow;
 617       borrow = y >> 16;
 618       Sign_Extend (borrow, y);
 619       *xc++ = y & 0xffff;
 620     }
 621   while (xb < xbe);
 622   while (xa < xae)
 623     {
 624       y = *xa++ + borrow;
 625       borrow = y >> 16;
 626       Sign_Extend (borrow, y);
 627       *xc++ = y & 0xffff;
 628     }
 629 #endif
 630   while (!*--xc)
 631     wa--;
 632   c->_wds = wa;
 633   return c;
 634 }
 635
 636 double
 637 ulp (double _x)
 638 {
 639   union double_union x, a;
 640   register __Long L;
 641
 642   x.d = _x;
 643
 644   L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
 645 #ifndef Sudden_Underflow
 646   if (L > 0)
 647     {
 648 #endif
 649 #ifdef IBM
 650       L |= Exp_msk1 >> 4;
 651 #endif
 652       word0 (a) = L;
 653 #ifndef _DOUBLE_IS_32BITS
 654       word1 (a) = 0;
 655 #endif
 656
 657 #ifndef Sudden_Underflow
 658     }
 659   else
 660     {
 661       L = -L >> Exp_shift;
 662       if (L < Exp_shift)
 663         {
 664           word0 (a) = 0x80000 >> L;
 665 #ifndef _DOUBLE_IS_32BITS
 666           word1 (a) = 0;
 667 #endif
 668         }
 669       else
 670         {
 671           word0 (a) = 0;
 672           L -= Exp_shift;
 673 #ifndef _DOUBLE_IS_32BITS
 674          word1 (a) = L >= 31 ? 1 : 1 << (31 - L);
 675 #endif
 676         }
 677     }
 678 #endif
 679   return a.d;
 680 }
 681
 682 double
 683 b2d (_Bigint * a, 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 d2b (struct _reent * ptr,
 760         double _d,
 761         int *e,
 762         int *bits)
 763
 764 {
 765   union double_union d;
 766   _Bigint *b;
 767   int de, i, k;
 768   __ULong *x, y, z;
 769 #ifdef VAX
 770   __ULong d0, d1;
 771 #endif
 772   d.d = _d;
 773 #ifdef VAX
 774   d0 = word0 (d) >> 16 | word0 (d) << 16;
 775   d1 = word1 (d) >> 16 | word1 (d) << 16;
 776 #else
 777 #define d0 word0(d)
 778 #define d1 word1(d)
 779   d.d = _d;
 780 #endif
 781
 782 #ifdef Pack_32
 783   b = eBalloc (ptr, 1);
 784 #else
 785   b = eBalloc (ptr, 2);
 786 #endif
 787   x = b->_x;
 788
 789   z = d0 & Frac_mask;
 790   d0 &= 0x7fffffff;             /* clear sign bit, which we ignore */
 791 #ifdef Sudden_Underflow
 792   de = (int) (d0 >> Exp_shift);
 793 #ifndef IBM
 794   z |= Exp_msk11;
 795 #endif
 796 #else
 797   if ((de = (int) (d0 >> Exp_shift)) != 0)
 798     z |= Exp_msk1;
 799 #endif
 800 #ifdef Pack_32
 801 #ifndef _DOUBLE_IS_32BITS
 802   if (d1)
 803     {
 804       y = d1;
 805       k = lo0bits (&y);
 806       if (k)
 807         {
 808          x[0] = y | z << (32 - k);
 809           z >>= k;
 810         }
 811       else
 812         x[0] = y;
 813       i = b->_wds = (x[1] = z) ? 2 : 1;
 814     }
 815   else
 816 #endif
 817     {
 818 #ifdef DEBUG
 819       if (!z)
 820         Bug ("Zero passed to d2b");
 821 #endif
 822       k = lo0bits (&z);
 823       x[0] = z;
 824       i = b->_wds = 1;
 825 #ifndef _DOUBLE_IS_32BITS
 826       k += 32;
 827 #endif
 828     }
 829 #else
 830   if (d1)
 831     {
 832       y = d1;
 833       k = lo0bits (&y);
 834       if (k)
 835         if (k >= 16)
 836           {
 837             x[0] = y | z << 32 - k & 0xffff;
 838             x[1] = z >> k - 16 & 0xffff;
 839             x[2] = z >> k;
 840             i = 2;
 841           }
 842         else
 843           {
 844             x[0] = y & 0xffff;
 845             x[1] = y >> 16 | z << 16 - k & 0xffff;
 846             x[2] = z >> k & 0xffff;
 847             x[3] = z >> k + 16;
 848             i = 3;
 849           }
 850       else
 851         {
 852           x[0] = y & 0xffff;
 853           x[1] = y >> 16;
 854           x[2] = z & 0xffff;
 855           x[3] = z >> 16;
 856           i = 3;
 857         }
 858     }
 859   else
 860     {
 861 #ifdef DEBUG
 862       if (!z)
 863         Bug ("Zero passed to d2b");
 864 #endif
 865       k = lo0bits (&z);
 866       if (k >= 16)
 867         {
 868           x[0] = z;
 869           i = 0;
 870         }
 871       else
 872         {
 873           x[0] = z & 0xffff;
 874           x[1] = z >> 16;
 875           i = 1;
 876         }
 877       k += 32;
 878     }
 879   while (!x[i])
 880     --i;
 881   b->_wds = i + 1;
 882 #endif
 883 #ifndef Sudden_Underflow
 884   if (de)
 885     {
 886 #endif
 887 #ifdef IBM
 888       *e = (de - Bias - (P - 1) << 2) + k;
 889       *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
 890 #else
 891       *e = de - Bias - (P - 1) + k;
 892       *bits = P - k;
 893 #endif
 894 #ifndef Sudden_Underflow
 895     }
 896   else
 897     {
 898       *e = de - Bias - (P - 1) + 1 + k;
 899 #ifdef Pack_32
 900       *bits = 32 * i - hi0bits (x[i - 1]);
 901 #else
 902       *bits = (i + 2) * 16 - hi0bits (x[i]);
 903 #endif
 904     }
 905 #endif
 906   return b;
 907 }
 908 #undef d0
 909 #undef d1
 910
 911 double
 912 ratio (_Bigint * a, _Bigint * b)
 913
 914 {
 915   union double_union da, db;
 916   int k, ka, kb;
 917
 918   da.d = b2d (a, &ka);
 919   db.d = b2d (b, &kb);
 920 #ifdef Pack_32
 921   k = ka - kb + 32 * (a->_wds - b->_wds);
 922 #else
 923   k = ka - kb + 16 * (a->_wds - b->_wds);
 924 #endif
 925 #ifdef IBM
 926   if (k > 0)
 927     {
 928       word0 (da) += (k >> 2) * Exp_msk1;
 929       if (k &= 3)
 930         da.d *= 1 << k;
 931     }
 932   else
 933     {
 934       k = -k;
 935       word0 (db) += (k >> 2) * Exp_msk1;
 936       if (k &= 3)
 937         db.d *= 1 << k;
 938     }
 939 #else
 940   if (k > 0)
 941     word0 (da) += k * Exp_msk1;
 942   else
 943     {
 944       k = -k;
 945       word0 (db) += k * Exp_msk1;
 946     }
 947 #endif
 948   return da.d / db.d;
 949 }
 950
 951
 952 const double
 953   tens[] =
 954 {
 955   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
 956   1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
 957   1e20, 1e21, 1e22, 1e23, 1e24
 958
 959 };
 960
 961 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
 962 const double bigtens[] =
 963 {1e16, 1e32, 1e64, 1e128, 1e256};
 964
 965 const double tinytens[] =
 966 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
 967 #else
 968 const double bigtens[] =
 969 {1e16, 1e32};
 970
 971 const double tinytens[] =
 972 {1e-16, 1e-32};
 973 #endif
 974
 975
 976 double
 977 _mprec_log10 (int dig)
 978 {
 979   double v = 1.0;
 980   if (dig < 24)
 981     return tens[dig];
 982   while (dig > 0)
 983     {
 984       v *= 10;
 985       dig--;
 986     }
 987   return v;
 988 }
 989
 990 void
 991 copybits (__ULong *c,
 992         int n,
 993         _Bigint *b)
 994 {
 995         __ULong *ce, *x, *xe;
 996 #ifdef Pack_16
 997         int nw, nw1;
 998 #endif
 999
1000         ce = c + ((n-1) >> kshift) + 1;
1001         x = b->_x;
1002 #ifdef Pack_32
1003         xe = x + b->_wds;
1004         while(x < xe)
1005                 *c++ = *x++;
1006 #else
1007         nw = b->_wds;
1008         nw1 = nw & 1;
1009         for(xe = x + (nw - nw1); x < xe; x += 2)
1010                 Storeinc(c, x[1], x[0]);
1011         if (nw1)
1012                 *c++ = *x;
1013 #endif
1014         while(c < ce)
1015                 *c++ = 0;
1016 }
1017
1018 __ULong
1019 any_on (_Bigint *b,
1020         int k)
1021 {
1022         int n, nwds;
1023         __ULong *x, *x0, x1, x2;
1024
1025         x = b->_x;
1026         nwds = b->_wds;
1027         n = k >> kshift;
1028         if (n > nwds)
1029                 n = nwds;
1030         else if (n < nwds && (k &= kmask)) {
1031                 x1 = x2 = x[n];
1032                 x1 >>= k;
1033                 x1 <<= k;
1034                 if (x1 != x2)
1035                         return 1;
1036                 }
1037         x0 = x;
1038         x += n;
1039         while(x > x0)
1040                 if (*--x)
1041                         return 1;
1042         return 0;
1043 }
1044