| /**************************************************************** |
| * |
| * The author of this software is David M. Gay. |
| * |
| * Copyright (c) 1991 by AT&T. |
| * |
| * Permission to use, copy, modify, and distribute this software for any |
| * purpose without fee is hereby granted, provided that this entire notice |
| * is included in all copies of any software which is or includes a copy |
| * or modification of this software and in all copies of the supporting |
| * documentation for such software. |
| * |
| * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |
| * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY |
| * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY |
| * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |
| * |
| ***************************************************************/ |
| |
| /* Please send bug reports to |
| David M. Gay |
| AT&T Bell Laboratories, Room 2C-463 |
| 600 Mountain Avenue |
| Murray Hill, NJ 07974-2070 |
| U.S.A. |
| dmg@research.att.com or research!dmg |
| */ |
| |
| /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. |
| * |
| * This strtod returns a nearest machine number to the input decimal |
| * string (or sets errno to ERANGE). With IEEE arithmetic, ties are |
| * broken by the IEEE round-even rule. Otherwise ties are broken by |
| * biased rounding (add half and chop). |
| * |
| * Inspired loosely by William D. Clinger's paper "How to Read Floating |
| * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. |
| * |
| * Modifications: |
| * |
| * 1. We only require IEEE, IBM, or VAX double-precision |
| * arithmetic (not IEEE double-extended). |
| * 2. We get by with floating-point arithmetic in a case that |
| * Clinger missed -- when we're computing d * 10^n |
| * for a small integer d and the integer n is not too |
| * much larger than 22 (the maximum integer k for which |
| * we can represent 10^k exactly), we may be able to |
| * compute (d*10^k) * 10^(e-k) with just one roundoff. |
| * 3. Rather than a bit-at-a-time adjustment of the binary |
| * result in the hard case, we use floating-point |
| * arithmetic to determine the adjustment to within |
| * one bit; only in really hard cases do we need to |
| * compute a second residual. |
| * 4. Because of 3., we don't need a large table of powers of 10 |
| * for ten-to-e (just some small tables, e.g. of 10^k |
| * for 0 <= k <= 22). |
| */ |
| |
| /* |
| * #define IEEE_8087 for IEEE-arithmetic machines where the least |
| * significant byte has the lowest address. |
| * #define IEEE_MC68k for IEEE-arithmetic machines where the most |
| * significant byte has the lowest address. |
| * #define Sudden_Underflow for IEEE-format machines without gradual |
| * underflow (i.e., that flush to zero on underflow). |
| * #define IBM for IBM mainframe-style floating-point arithmetic. |
| * #define VAX for VAX-style floating-point arithmetic. |
| * #define Unsigned_Shifts if >> does treats its left operand as unsigned. |
| * #define No_leftright to omit left-right logic in fast floating-point |
| * computation of dtoa. |
| * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. |
| * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines |
| * that use extended-precision instructions to compute rounded |
| * products and quotients) with IBM. |
| * #define ROUND_BIASED for IEEE-format with biased rounding. |
| * #define Inaccurate_Divide for IEEE-format with correctly rounded |
| * products but inaccurate quotients, e.g., for Intel i860. |
| * #define Just_16 to store 16 bits per 32-bit long when doing high-precision |
| * integer arithmetic. Whether this speeds things up or slows things |
| * down depends on the machine and the number being converted. |
| */ |
| |
| #include <stdlib.h> |
| #include <string.h> |
| #include <java-assert.h> |
| #include "mprec.h" |
| |
| /* reent.c knows this value */ |
| #define _Kmax 15 |
| #include <stdio.h> |
| |
| _Jv_Bigint * |
| _DEFUN (Balloc, (ptr, k), struct _Jv_reent *ptr _AND int k) |
| { |
| _Jv_Bigint *rv = NULL; |
| |
| int i = 0; |
| int j = 1; |
| |
| JvAssert ((1 << k) < MAX_BIGNUM_WDS); |
| |
| while ((ptr->_allocation_map & j) && i < MAX_BIGNUMS) |
| i++, j <<= 1; |
| |
| JvAssert (i < MAX_BIGNUMS); |
| |
| if (i >= MAX_BIGNUMS) |
| return NULL; |
| |
| ptr->_allocation_map |= j; |
| rv = &ptr->_freelist[i]; |
| |
| rv->_k = k; |
| rv->_maxwds = 32; |
| |
| return rv; |
| } |
| |
| |
| void |
| _DEFUN (Bfree, (ptr, v), struct _Jv_reent *ptr _AND _Jv_Bigint * v) |
| { |
| long i; |
| |
| i = v - ptr->_freelist; |
| |
| JvAssert (i >= 0 && i < MAX_BIGNUMS); |
| |
| if (i >= 0 && i < MAX_BIGNUMS) |
| ptr->_allocation_map &= ~ (1 << i); |
| } |
| |
| |
| _Jv_Bigint * |
| _DEFUN (multadd, (ptr, b, m, a), |
| struct _Jv_reent *ptr _AND |
| _Jv_Bigint * b _AND |
| int m _AND |
| int a) |
| { |
| int i, wds; |
| unsigned long *x, y; |
| #ifdef Pack_32 |
| unsigned long xi, z; |
| #endif |
| _Jv_Bigint *b1; |
| |
| wds = b->_wds; |
| x = b->_x; |
| i = 0; |
| do |
| { |
| #ifdef Pack_32 |
| xi = *x; |
| y = (xi & 0xffff) * m + a; |
| z = (xi >> 16) * m + (y >> 16); |
| a = (int) (z >> 16); |
| *x++ = (z << 16) + (y & 0xffff); |
| #else |
| y = *x * m + a; |
| a = (int) (y >> 16); |
| *x++ = y & 0xffff; |
| #endif |
| } |
| while (++i < wds); |
| if (a) |
| { |
| if (wds >= b->_maxwds) |
| { |
| b1 = Balloc (ptr, b->_k + 1); |
| Bcopy (b1, b); |
| Bfree (ptr, b); |
| b = b1; |
| } |
| b->_x[wds++] = a; |
| b->_wds = wds; |
| } |
| return b; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (s2b, (ptr, s, nd0, nd, y9), |
| struct _Jv_reent * ptr _AND |
| _CONST char *s _AND |
| int nd0 _AND |
| int nd _AND |
| unsigned long y9) |
| { |
| _Jv_Bigint *b; |
| int i, k; |
| long x, y; |
| |
| x = (nd + 8) / 9; |
| for (k = 0, y = 1; x > y; y <<= 1, k++); |
| #ifdef Pack_32 |
| b = Balloc (ptr, k); |
| b->_x[0] = y9; |
| b->_wds = 1; |
| #else |
| b = Balloc (ptr, k + 1); |
| b->_x[0] = y9 & 0xffff; |
| b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1; |
| #endif |
| |
| i = 9; |
| if (9 < nd0) |
| { |
| s += 9; |
| do |
| b = multadd (ptr, b, 10, *s++ - '0'); |
| while (++i < nd0); |
| s++; |
| } |
| else |
| s += 10; |
| for (; i < nd; i++) |
| b = multadd (ptr, b, 10, *s++ - '0'); |
| return b; |
| } |
| |
| int |
| _DEFUN (hi0bits, |
| (x), register unsigned long x) |
| { |
| register int k = 0; |
| |
| if (!(x & 0xffff0000)) |
| { |
| k = 16; |
| x <<= 16; |
| } |
| if (!(x & 0xff000000)) |
| { |
| k += 8; |
| x <<= 8; |
| } |
| if (!(x & 0xf0000000)) |
| { |
| k += 4; |
| x <<= 4; |
| } |
| if (!(x & 0xc0000000)) |
| { |
| k += 2; |
| x <<= 2; |
| } |
| if (!(x & 0x80000000)) |
| { |
| k++; |
| if (!(x & 0x40000000)) |
| return 32; |
| } |
| return k; |
| } |
| |
| int |
| _DEFUN (lo0bits, (y), unsigned long *y) |
| { |
| register int k; |
| register unsigned long x = *y; |
| |
| if (x & 7) |
| { |
| if (x & 1) |
| return 0; |
| if (x & 2) |
| { |
| *y = x >> 1; |
| return 1; |
| } |
| *y = x >> 2; |
| return 2; |
| } |
| k = 0; |
| if (!(x & 0xffff)) |
| { |
| k = 16; |
| x >>= 16; |
| } |
| if (!(x & 0xff)) |
| { |
| k += 8; |
| x >>= 8; |
| } |
| if (!(x & 0xf)) |
| { |
| k += 4; |
| x >>= 4; |
| } |
| if (!(x & 0x3)) |
| { |
| k += 2; |
| x >>= 2; |
| } |
| if (!(x & 1)) |
| { |
| k++; |
| x >>= 1; |
| if (!(x & 1)) |
| return 32; |
| } |
| *y = x; |
| return k; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (i2b, (ptr, i), struct _Jv_reent * ptr _AND int i) |
| { |
| _Jv_Bigint *b; |
| |
| b = Balloc (ptr, 1); |
| b->_x[0] = i; |
| b->_wds = 1; |
| return b; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (mult, (ptr, a, b), struct _Jv_reent * ptr _AND _Jv_Bigint * a _AND _Jv_Bigint * b) |
| { |
| _Jv_Bigint *c; |
| int k, wa, wb, wc; |
| unsigned long carry, y, z; |
| unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0; |
| #ifdef Pack_32 |
| unsigned long z2; |
| #endif |
| |
| if (a->_wds < b->_wds) |
| { |
| c = a; |
| a = b; |
| b = c; |
| } |
| k = a->_k; |
| wa = a->_wds; |
| wb = b->_wds; |
| wc = wa + wb; |
| if (wc > a->_maxwds) |
| k++; |
| c = Balloc (ptr, k); |
| for (x = c->_x, xa = x + wc; x < xa; x++) |
| *x = 0; |
| xa = a->_x; |
| xae = xa + wa; |
| xb = b->_x; |
| xbe = xb + wb; |
| xc0 = c->_x; |
| #ifdef Pack_32 |
| for (; xb < xbe; xb++, xc0++) |
| { |
| if ((y = *xb & 0xffff)) |
| { |
| x = xa; |
| xc = xc0; |
| carry = 0; |
| do |
| { |
| z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; |
| carry = z >> 16; |
| z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; |
| carry = z2 >> 16; |
| Storeinc (xc, z2, z); |
| } |
| while (x < xae); |
| *xc = carry; |
| } |
| if ((y = *xb >> 16)) |
| { |
| x = xa; |
| xc = xc0; |
| carry = 0; |
| z2 = *xc; |
| do |
| { |
| z = (*x & 0xffff) * y + (*xc >> 16) + carry; |
| carry = z >> 16; |
| Storeinc (xc, z, z2); |
| z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; |
| carry = z2 >> 16; |
| } |
| while (x < xae); |
| *xc = z2; |
| } |
| } |
| #else |
| for (; xb < xbe; xc0++) |
| { |
| if (y = *xb++) |
| { |
| x = xa; |
| xc = xc0; |
| carry = 0; |
| do |
| { |
| z = *x++ * y + *xc + carry; |
| carry = z >> 16; |
| *xc++ = z & 0xffff; |
| } |
| while (x < xae); |
| *xc = carry; |
| } |
| } |
| #endif |
| for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc); |
| c->_wds = wc; |
| return c; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (pow5mult, |
| (ptr, b, k), struct _Jv_reent * ptr _AND _Jv_Bigint * b _AND int k) |
| { |
| _Jv_Bigint *b1, *p5, *p51; |
| int i; |
| static _CONST int p05[3] = {5, 25, 125}; |
| |
| if ((i = k & 3)) |
| b = multadd (ptr, b, p05[i - 1], 0); |
| |
| if (!(k >>= 2)) |
| return b; |
| if (!(p5 = ptr->_p5s)) |
| { |
| /* first time */ |
| p5 = ptr->_p5s = i2b (ptr, 625); |
| p5->_next = 0; |
| } |
| for (;;) |
| { |
| if (k & 1) |
| { |
| b1 = mult (ptr, b, p5); |
| Bfree (ptr, b); |
| b = b1; |
| } |
| if (!(k >>= 1)) |
| break; |
| if (!(p51 = p5->_next)) |
| { |
| p51 = p5->_next = mult (ptr, p5, p5); |
| p51->_next = 0; |
| } |
| p5 = p51; |
| } |
| return b; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (lshift, (ptr, b, k), struct _Jv_reent * ptr _AND _Jv_Bigint * b _AND int k) |
| { |
| int i, k1, n, n1; |
| _Jv_Bigint *b1; |
| unsigned long *x, *x1, *xe, z; |
| |
| #ifdef Pack_32 |
| n = k >> 5; |
| #else |
| n = k >> 4; |
| #endif |
| k1 = b->_k; |
| n1 = n + b->_wds + 1; |
| for (i = b->_maxwds; n1 > i; i <<= 1) |
| k1++; |
| b1 = Balloc (ptr, k1); |
| x1 = b1->_x; |
| for (i = 0; i < n; i++) |
| *x1++ = 0; |
| x = b->_x; |
| xe = x + b->_wds; |
| #ifdef Pack_32 |
| if (k &= 0x1f) |
| { |
| k1 = 32 - k; |
| z = 0; |
| do |
| { |
| *x1++ = *x << k | z; |
| z = *x++ >> k1; |
| } |
| while (x < xe); |
| if ((*x1 = z)) |
| ++n1; |
| } |
| #else |
| if (k &= 0xf) |
| { |
| k1 = 16 - k; |
| z = 0; |
| do |
| { |
| *x1++ = *x << k & 0xffff | z; |
| z = *x++ >> k1; |
| } |
| while (x < xe); |
| if (*x1 = z) |
| ++n1; |
| } |
| #endif |
| else |
| do |
| *x1++ = *x++; |
| while (x < xe); |
| b1->_wds = n1 - 1; |
| Bfree (ptr, b); |
| return b1; |
| } |
| |
| int |
| _DEFUN (cmp, (a, b), _Jv_Bigint * a _AND _Jv_Bigint * b) |
| { |
| unsigned long *xa, *xa0, *xb, *xb0; |
| int i, j; |
| |
| i = a->_wds; |
| j = b->_wds; |
| #ifdef DEBUG |
| if (i > 1 && !a->_x[i - 1]) |
| Bug ("cmp called with a->_x[a->_wds-1] == 0"); |
| if (j > 1 && !b->_x[j - 1]) |
| Bug ("cmp called with b->_x[b->_wds-1] == 0"); |
| #endif |
| if (i -= j) |
| return i; |
| xa0 = a->_x; |
| xa = xa0 + j; |
| xb0 = b->_x; |
| xb = xb0 + j; |
| for (;;) |
| { |
| if (*--xa != *--xb) |
| return *xa < *xb ? -1 : 1; |
| if (xa <= xa0) |
| break; |
| } |
| return 0; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (diff, (ptr, a, b), struct _Jv_reent * ptr _AND |
| _Jv_Bigint * a _AND _Jv_Bigint * b) |
| { |
| _Jv_Bigint *c; |
| int i, wa, wb; |
| long borrow, y; /* We need signed shifts here. */ |
| unsigned long *xa, *xae, *xb, *xbe, *xc; |
| #ifdef Pack_32 |
| long z; |
| #endif |
| |
| i = cmp (a, b); |
| if (!i) |
| { |
| c = Balloc (ptr, 0); |
| c->_wds = 1; |
| c->_x[0] = 0; |
| return c; |
| } |
| if (i < 0) |
| { |
| c = a; |
| a = b; |
| b = c; |
| i = 1; |
| } |
| else |
| i = 0; |
| c = Balloc (ptr, a->_k); |
| c->_sign = i; |
| wa = a->_wds; |
| xa = a->_x; |
| xae = xa + wa; |
| wb = b->_wds; |
| xb = b->_x; |
| xbe = xb + wb; |
| xc = c->_x; |
| borrow = 0; |
| #ifdef Pack_32 |
| do |
| { |
| y = (*xa & 0xffff) - (*xb & 0xffff) + borrow; |
| borrow = y >> 16; |
| Sign_Extend (borrow, y); |
| z = (*xa++ >> 16) - (*xb++ >> 16) + borrow; |
| borrow = z >> 16; |
| Sign_Extend (borrow, z); |
| Storeinc (xc, z, y); |
| } |
| while (xb < xbe); |
| while (xa < xae) |
| { |
| y = (*xa & 0xffff) + borrow; |
| borrow = y >> 16; |
| Sign_Extend (borrow, y); |
| z = (*xa++ >> 16) + borrow; |
| borrow = z >> 16; |
| Sign_Extend (borrow, z); |
| Storeinc (xc, z, y); |
| } |
| #else |
| do |
| { |
| y = *xa++ - *xb++ + borrow; |
| borrow = y >> 16; |
| Sign_Extend (borrow, y); |
| *xc++ = y & 0xffff; |
| } |
| while (xb < xbe); |
| while (xa < xae) |
| { |
| y = *xa++ + borrow; |
| borrow = y >> 16; |
| Sign_Extend (borrow, y); |
| *xc++ = y & 0xffff; |
| } |
| #endif |
| while (!*--xc) |
| wa--; |
| c->_wds = wa; |
| return c; |
| } |
| |
| double |
| _DEFUN (ulp, (_x), double _x) |
| { |
| union double_union x, a; |
| register long L; |
| |
| x.d = _x; |
| |
| L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1; |
| #ifndef Sudden_Underflow |
| if (L > 0) |
| { |
| #endif |
| #ifdef IBM |
| L |= Exp_msk1 >> 4; |
| #endif |
| word0 (a) = L; |
| #ifndef _DOUBLE_IS_32BITS |
| word1 (a) = 0; |
| #endif |
| |
| #ifndef Sudden_Underflow |
| } |
| else |
| { |
| L = -L >> Exp_shift; |
| if (L < Exp_shift) |
| { |
| word0 (a) = 0x80000 >> L; |
| #ifndef _DOUBLE_IS_32BITS |
| word1 (a) = 0; |
| #endif |
| } |
| else |
| { |
| word0 (a) = 0; |
| L -= Exp_shift; |
| #ifndef _DOUBLE_IS_32BITS |
| word1 (a) = L >= 31 ? 1 : 1 << (31 - L); |
| #endif |
| } |
| } |
| #endif |
| return a.d; |
| } |
| |
| double |
| _DEFUN (b2d, (a, e), |
| _Jv_Bigint * a _AND int *e) |
| { |
| unsigned long *xa, *xa0, w, y, z; |
| int k; |
| union double_union d; |
| #ifdef VAX |
| unsigned long d0, d1; |
| #else |
| #define d0 word0(d) |
| #define d1 word1(d) |
| #endif |
| |
| xa0 = a->_x; |
| xa = xa0 + a->_wds; |
| y = *--xa; |
| #ifdef DEBUG |
| if (!y) |
| Bug ("zero y in b2d"); |
| #endif |
| k = hi0bits (y); |
| *e = 32 - k; |
| #ifdef Pack_32 |
| if (k < Ebits) |
| { |
| d0 = Exp_1 | y >> (Ebits - k); |
| w = xa > xa0 ? *--xa : 0; |
| #ifndef _DOUBLE_IS_32BITS |
| d1 = y << (32 - Ebits + k) | w >> (Ebits - k); |
| #endif |
| goto ret_d; |
| } |
| z = xa > xa0 ? *--xa : 0; |
| if (k -= Ebits) |
| { |
| d0 = Exp_1 | y << k | z >> (32 - k); |
| y = xa > xa0 ? *--xa : 0; |
| #ifndef _DOUBLE_IS_32BITS |
| d1 = z << k | y >> (32 - k); |
| #endif |
| } |
| else |
| { |
| d0 = Exp_1 | y; |
| #ifndef _DOUBLE_IS_32BITS |
| d1 = z; |
| #endif |
| } |
| #else |
| if (k < Ebits + 16) |
| { |
| z = xa > xa0 ? *--xa : 0; |
| d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; |
| w = xa > xa0 ? *--xa : 0; |
| y = xa > xa0 ? *--xa : 0; |
| d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; |
| goto ret_d; |
| } |
| z = xa > xa0 ? *--xa : 0; |
| w = xa > xa0 ? *--xa : 0; |
| k -= Ebits + 16; |
| d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; |
| y = xa > xa0 ? *--xa : 0; |
| d1 = w << k + 16 | y << k; |
| #endif |
| ret_d: |
| #ifdef VAX |
| word0 (d) = d0 >> 16 | d0 << 16; |
| word1 (d) = d1 >> 16 | d1 << 16; |
| #else |
| #undef d0 |
| #undef d1 |
| #endif |
| return d.d; |
| } |
| |
| _Jv_Bigint * |
| _DEFUN (d2b, |
| (ptr, _d, e, bits), |
| struct _Jv_reent * ptr _AND |
| double _d _AND |
| int *e _AND |
| int *bits) |
| |
| { |
| union double_union d; |
| _Jv_Bigint *b; |
| int de, i, k; |
| unsigned long *x, y, z; |
| #ifdef VAX |
| unsigned long d0, d1; |
| d.d = _d; |
| d0 = word0 (d) >> 16 | word0 (d) << 16; |
| d1 = word1 (d) >> 16 | word1 (d) << 16; |
| #else |
| #define d0 word0(d) |
| #define d1 word1(d) |
| d.d = _d; |
| #endif |
| |
| #ifdef Pack_32 |
| b = Balloc (ptr, 1); |
| #else |
| b = Balloc (ptr, 2); |
| #endif |
| x = b->_x; |
| |
| z = d0 & Frac_mask; |
| d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ |
| #ifdef Sudden_Underflow |
| de = (int) (d0 >> Exp_shift); |
| #ifndef IBM |
| z |= Exp_msk11; |
| #endif |
| #else |
| if ((de = (int) (d0 >> Exp_shift))) |
| z |= Exp_msk1; |
| #endif |
| #ifdef Pack_32 |
| #ifndef _DOUBLE_IS_32BITS |
| if ((y = d1)) |
| { |
| if ((k = lo0bits (&y))) |
| { |
| x[0] = y | z << (32 - k); |
| z >>= k; |
| } |
| else |
| x[0] = y; |
| i = b->_wds = (x[1] = z) ? 2 : 1; |
| } |
| else |
| #endif |
| { |
| #ifdef DEBUG |
| if (!z) |
| Bug ("Zero passed to d2b"); |
| #endif |
| k = lo0bits (&z); |
| x[0] = z; |
| i = b->_wds = 1; |
| #ifndef _DOUBLE_IS_32BITS |
| k += 32; |
| #endif |
| } |
| #else |
| if (y = d1) |
| { |
| if (k = lo0bits (&y)) |
| if (k >= 16) |
| { |
| x[0] = y | z << 32 - k & 0xffff; |
| x[1] = z >> k - 16 & 0xffff; |
| x[2] = z >> k; |
| i = 2; |
| } |
| else |
| { |
| x[0] = y & 0xffff; |
| x[1] = y >> 16 | z << 16 - k & 0xffff; |
| x[2] = z >> k & 0xffff; |
| x[3] = z >> k + 16; |
| i = 3; |
| } |
| else |
| { |
| x[0] = y & 0xffff; |
| x[1] = y >> 16; |
| x[2] = z & 0xffff; |
| x[3] = z >> 16; |
| i = 3; |
| } |
| } |
| else |
| { |
| #ifdef DEBUG |
| if (!z) |
| Bug ("Zero passed to d2b"); |
| #endif |
| k = lo0bits (&z); |
| if (k >= 16) |
| { |
| x[0] = z; |
| i = 0; |
| } |
| else |
| { |
| x[0] = z & 0xffff; |
| x[1] = z >> 16; |
| i = 1; |
| } |
| k += 32; |
| } |
| while (!x[i]) |
| --i; |
| b->_wds = i + 1; |
| #endif |
| #ifndef Sudden_Underflow |
| if (de) |
| { |
| #endif |
| #ifdef IBM |
| *e = (de - Bias - (P - 1) << 2) + k; |
| *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask); |
| #else |
| *e = de - Bias - (P - 1) + k; |
| *bits = P - k; |
| #endif |
| #ifndef Sudden_Underflow |
| } |
| else |
| { |
| *e = de - Bias - (P - 1) + 1 + k; |
| #ifdef Pack_32 |
| *bits = 32 * i - hi0bits (x[i - 1]); |
| #else |
| *bits = (i + 2) * 16 - hi0bits (x[i]); |
| #endif |
| } |
| #endif |
| return b; |
| } |
| #undef d0 |
| #undef d1 |
| |
| double |
| _DEFUN (ratio, (a, b), _Jv_Bigint * a _AND _Jv_Bigint * b) |
| |
| { |
| union double_union da, db; |
| int k, ka, kb; |
| |
| da.d = b2d (a, &ka); |
| db.d = b2d (b, &kb); |
| #ifdef Pack_32 |
| k = ka - kb + 32 * (a->_wds - b->_wds); |
| #else |
| k = ka - kb + 16 * (a->_wds - b->_wds); |
| #endif |
| #ifdef IBM |
| if (k > 0) |
| { |
| word0 (da) += (k >> 2) * Exp_msk1; |
| if (k &= 3) |
| da.d *= 1 << k; |
| } |
| else |
| { |
| k = -k; |
| word0 (db) += (k >> 2) * Exp_msk1; |
| if (k &= 3) |
| db.d *= 1 << k; |
| } |
| #else |
| if (k > 0) |
| word0 (da) += k * Exp_msk1; |
| else |
| { |
| k = -k; |
| word0 (db) += k * Exp_msk1; |
| } |
| #endif |
| return da.d / db.d; |
| } |
| |
| |
| _CONST double |
| tens[] = |
| { |
| 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, |
| 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, |
| 1e20, 1e21, 1e22, 1e23, 1e24 |
| |
| }; |
| |
| #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800) |
| _CONST double bigtens[] = |
| {1e16, 1e32, 1e64, 1e128, 1e256}; |
| |
| _CONST double tinytens[] = |
| {1e-16, 1e-32, 1e-64, 1e-128, 1e-256}; |
| #else |
| _CONST double bigtens[] = |
| {1e16, 1e32}; |
| |
| _CONST double tinytens[] = |
| {1e-16, 1e-32}; |
| #endif |
| |
| |