0,0 → 1,1049 |
/**************************************************************** |
* |
* 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 <_ansi.h> |
#include <stdlib.h> |
#include <string.h> |
#include <reent.h> |
#include "mprec.h" |
|
/* This is defined in sys/reent.h as (sizeof (size_t) << 3) now, as in NetBSD. |
The old value of 15 was wrong and made newlib vulnerable against buffer |
overrun attacks (CVE-2009-0689), same as other implementations of gdtoa |
based on BSD code. |
#define _Kmax 15 |
*/ |
|
_Bigint * |
_DEFUN (Balloc, (ptr, k), struct _reent *ptr _AND int k) |
{ |
int x; |
_Bigint *rv ; |
|
_REENT_CHECK_MP(ptr); |
if (_REENT_MP_FREELIST(ptr) == NULL) |
{ |
/* Allocate a list of pointers to the mprec objects */ |
_REENT_MP_FREELIST(ptr) = (struct _Bigint **) _calloc_r (ptr, |
sizeof (struct _Bigint *), |
_Kmax + 1); |
if (_REENT_MP_FREELIST(ptr) == NULL) |
{ |
return NULL; |
} |
} |
|
if ((rv = _REENT_MP_FREELIST(ptr)[k]) != 0) |
{ |
_REENT_MP_FREELIST(ptr)[k] = rv->_next; |
} |
else |
{ |
x = 1 << k; |
/* Allocate an mprec Bigint and stick in in the freelist */ |
rv = (_Bigint *) _calloc_r (ptr, |
1, |
sizeof (_Bigint) + |
(x-1) * sizeof(rv->_x)); |
if (rv == NULL) return NULL; |
rv->_k = k; |
rv->_maxwds = x; |
} |
rv->_sign = rv->_wds = 0; |
return rv; |
} |
|
void |
_DEFUN (Bfree, (ptr, v), struct _reent *ptr _AND _Bigint * v) |
{ |
_REENT_CHECK_MP(ptr); |
if (v) |
{ |
v->_next = _REENT_MP_FREELIST(ptr)[v->_k]; |
_REENT_MP_FREELIST(ptr)[v->_k] = v; |
} |
} |
|
_Bigint * |
_DEFUN (multadd, (ptr, b, m, a), |
struct _reent *ptr _AND |
_Bigint * b _AND |
int m _AND |
int a) |
{ |
int i, wds; |
__ULong *x, y; |
#ifdef Pack_32 |
__ULong xi, z; |
#endif |
_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; |
} |
|
_Bigint * |
_DEFUN (s2b, (ptr, s, nd0, nd, y9), |
struct _reent * ptr _AND |
_CONST char *s _AND |
int nd0 _AND |
int nd _AND |
__ULong y9) |
{ |
_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 __ULong 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), __ULong *y) |
{ |
register int k; |
register __ULong 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; |
} |
|
_Bigint * |
_DEFUN (i2b, (ptr, i), struct _reent * ptr _AND int i) |
{ |
_Bigint *b; |
|
b = Balloc (ptr, 1); |
b->_x[0] = i; |
b->_wds = 1; |
return b; |
} |
|
_Bigint * |
_DEFUN (mult, (ptr, a, b), struct _reent * ptr _AND _Bigint * a _AND _Bigint * b) |
{ |
_Bigint *c; |
int k, wa, wb, wc; |
__ULong carry, y, z; |
__ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; |
#ifdef Pack_32 |
__ULong 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) != 0) |
{ |
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) != 0) |
{ |
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; |
} |
|
_Bigint * |
_DEFUN (pow5mult, |
(ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) |
{ |
_Bigint *b1, *p5, *p51; |
int i; |
static _CONST int p05[3] = {5, 25, 125}; |
|
if ((i = k & 3) != 0) |
b = multadd (ptr, b, p05[i - 1], 0); |
|
if (!(k >>= 2)) |
return b; |
_REENT_CHECK_MP(ptr); |
if (!(p5 = _REENT_MP_P5S(ptr))) |
{ |
/* first time */ |
p5 = _REENT_MP_P5S(ptr) = 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; |
} |
|
_Bigint * |
_DEFUN (lshift, (ptr, b, k), struct _reent * ptr _AND _Bigint * b _AND int k) |
{ |
int i, k1, n, n1; |
_Bigint *b1; |
__ULong *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) != 0) |
++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), _Bigint * a _AND _Bigint * b) |
{ |
__ULong *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; |
} |
|
_Bigint * |
_DEFUN (diff, (ptr, a, b), struct _reent * ptr _AND |
_Bigint * a _AND _Bigint * b) |
{ |
_Bigint *c; |
int i, wa, wb; |
__Long borrow, y; /* We need signed shifts here. */ |
__ULong *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), |
_Bigint * a _AND int *e) |
{ |
__ULong *xa, *xa0, w, y, z; |
int k; |
union double_union d; |
#ifdef VAX |
__ULong 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; |
} |
|
_Bigint * |
_DEFUN (d2b, |
(ptr, _d, e, bits), |
struct _reent * ptr _AND |
double _d _AND |
int *e _AND |
int *bits) |
|
{ |
union double_union d; |
_Bigint *b; |
int de, i, k; |
__ULong *x, y, z; |
#ifdef VAX |
__ULong d0, d1; |
#endif |
d.d = _d; |
#ifdef VAX |
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)) != 0) |
z |= Exp_msk1; |
#endif |
#ifdef Pack_32 |
#ifndef _DOUBLE_IS_32BITS |
if (d1) |
{ |
y = d1; |
k = lo0bits (&y); |
if (k) |
{ |
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 (d1) |
{ |
y = d1; |
k = lo0bits (&y); |
if (k) |
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), _Bigint * a _AND _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 |
|
|
double |
_DEFUN (_mprec_log10, (dig), |
int dig) |
{ |
double v = 1.0; |
if (dig < 24) |
return tens[dig]; |
while (dig > 0) |
{ |
v *= 10; |
dig--; |
} |
return v; |
} |
|
void |
_DEFUN (copybits, (c, n, b), |
__ULong *c _AND |
int n _AND |
_Bigint *b) |
{ |
__ULong *ce, *x, *xe; |
#ifdef Pack_16 |
int nw, nw1; |
#endif |
|
ce = c + ((n-1) >> kshift) + 1; |
x = b->_x; |
#ifdef Pack_32 |
xe = x + b->_wds; |
while(x < xe) |
*c++ = *x++; |
#else |
nw = b->_wds; |
nw1 = nw & 1; |
for(xe = x + (nw - nw1); x < xe; x += 2) |
Storeinc(c, x[1], x[0]); |
if (nw1) |
*c++ = *x; |
#endif |
while(c < ce) |
*c++ = 0; |
} |
|
__ULong |
_DEFUN (any_on, (b, k), |
_Bigint *b _AND |
int k) |
{ |
int n, nwds; |
__ULong *x, *x0, x1, x2; |
|
x = b->_x; |
nwds = b->_wds; |
n = k >> kshift; |
if (n > nwds) |
n = nwds; |
else if (n < nwds && (k &= kmask)) { |
x1 = x2 = x[n]; |
x1 >>= k; |
x1 <<= k; |
if (x1 != x2) |
return 1; |
} |
x0 = x; |
x += n; |
while(x > x0) |
if (*--x) |
return 1; |
return 0; |
} |
|