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