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