25 #ifdef WORDS_BIGENDIAN
36 #if SIZEOF_LONG == 2*SIZEOF_INT
38 #define Intcast (int)(long)
246 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
256 #ifdef Honor_FLT_ROUNDS
257 #ifndef Trust_FLT_ROUNDS
266 extern void *
MALLOC(
size_t);
269 #define MALLOC malloc
272 #ifndef Omit_Private_Memory
274 #define PRIVATE_MEM 2304
276 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
281 #undef Avoid_Underflow
290 #ifndef NO_INFNAN_CHECK
296 #define NO_STRTOD_BIGCOMP
305 #define DBL_MAX_10_EXP 308
306 #define DBL_MAX_EXP 1024
312 #define DBL_MAX_10_EXP 75
313 #define DBL_MAX_EXP 63
315 #define DBL_MAX 7.2370055773322621e+75
320 #define DBL_MAX_10_EXP 38
321 #define DBL_MAX_EXP 127
323 #define DBL_MAX 1.7014118346046923e+38
327 #define LONG_MAX 2147483647
357 #define word0(x) (x)->L[1]
358 #define word1(x) (x)->L[0]
360 #define word0(x) (x)->L[0]
361 #define word1(x) (x)->L[1]
363 #define dval(x) (x)->d
365 #ifndef STRTOD_DIGLIM
366 #define STRTOD_DIGLIM 40
372 #define strtod_diglim STRTOD_DIGLIM
379 #if defined(IEEE_8087) + defined(VAX)
380 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
381 ((unsigned short *)a)[0] = (unsigned short)c, a++)
383 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
384 ((unsigned short *)a)[1] = (unsigned short)c, a++)
395 #define Exp_shift1 20
396 #define Exp_msk1 0x100000
397 #define Exp_msk11 0x100000
398 #define Exp_mask 0x7ff00000
404 #define Exp_1 0x3ff00000
405 #define Exp_11 0x3ff00000
407 #define Frac_mask 0xfffff
408 #define Frac_mask1 0xfffff
411 #define Bndry_mask 0xfffff
412 #define Bndry_mask1 0xfffff
414 #define Sign_bit 0x80000000
420 #ifndef NO_IEEE_Scale
421 #define Avoid_Underflow
423 #undef Sudden_Underflow
429 #define Flt_Rounds FLT_ROUNDS
435 #ifdef Honor_FLT_ROUNDS
436 #undef Check_FLT_ROUNDS
437 #define Check_FLT_ROUNDS
439 #define Rounding Flt_Rounds
443 #undef Check_FLT_ROUNDS
444 #undef Honor_FLT_ROUNDS
446 #undef Sudden_Underflow
447 #define Sudden_Underflow
452 #define Exp_shift1 24
453 #define Exp_msk1 0x1000000
454 #define Exp_msk11 0x1000000
455 #define Exp_mask 0x7f000000
461 #define Exp_1 0x41000000
462 #define Exp_11 0x41000000
464 #define Frac_mask 0xffffff
465 #define Frac_mask1 0xffffff
468 #define Bndry_mask 0xefffff
469 #define Bndry_mask1 0xffffff
471 #define Sign_bit 0x80000000
473 #define Tiny0 0x100000
482 #define Exp_msk1 0x80
483 #define Exp_msk11 0x800000
484 #define Exp_mask 0x7f80
490 #define Exp_1 0x40800000
491 #define Exp_11 0x4080
493 #define Frac_mask 0x7fffff
494 #define Frac_mask1 0xffff007f
497 #define Bndry_mask 0xffff007f
498 #define Bndry_mask1 0xffff007f
500 #define Sign_bit 0x8000
512 #ifdef ROUND_BIASED_without_Round_Up
519 #define rounded_product(a,b) a = rnd_prod(a, b)
520 #define rounded_quotient(a,b) a = rnd_quot(a, b)
522 extern double rnd_prod(), rnd_quot();
524 extern double rnd_prod(
double,
double), rnd_quot(
double,
double);
527 #define rounded_product(a,b) a *= b
528 #define rounded_quotient(a,b) a /= b
531 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
532 #define Big1 0xffffffff
540 BCinfo {
int dp0,
dp1,
dplen,
dsign,
e0,
inexact,
nd,
nd0,
rounding,
scale,
uflchk; };
543 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
545 #define FFFFFFFF 0xffffffffUL
560 #define Llong long long
563 #define ULLong unsigned Llong
567 #ifndef MULTIPLE_THREADS
568 #define ACQUIRE_DTOA_LOCK(n)
569 #define FREE_DTOA_LOCK(n)
575 extern "C" double os_strtod(
const char *s00,
char **se);
576 extern "C" char *
os_dtoa(
double d,
int mode,
int ndigits,
577 int *decpt,
int *sign,
char **rve);
601 #ifndef Omit_Private_Memory
612 #ifdef Omit_Private_Memory
615 len = (
sizeof(
Bigint) + (
x-1)*
sizeof(
ULong) +
sizeof(
double) - 1)
656 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
657 y->wds*sizeof(Long) + 2*sizeof(int))
662 (b, m, a)
Bigint *b;
int m, a;
685 y = *
x * (ULLong)m + carry;
691 y = (xi & 0xffff) * m + carry;
692 z = (xi >> 16) * m + (y >> 16);
694 *
x++ = (z << 16) + (y & 0xffff);
704 if (
wds >= b->maxwds) {
719 (s, nd0, nd, y9, dplen)
CONST char *s;
int nd0, nd, dplen;
ULong y9;
721 (
const char *s,
int nd0,
int nd,
ULong y9,
int dplen)
729 for(
k = 0, y = 1;
x > y; y <<= 1,
k++) ;
736 b->
x[0] = y9 & 0xffff;
737 b->
wds = (b->
x[1] = y9 >> 16) ? 2 : 1;
743 do b =
multadd(b, 10, *s++ -
'0');
750 b =
multadd(b, 10, *s++ -
'0');
764 if (!(
x & 0xffff0000)) {
768 if (!(
x & 0xff000000)) {
772 if (!(
x & 0xf0000000)) {
776 if (!(
x & 0xc0000000)) {
780 if (!(
x & 0x80000000)) {
782 if (!(
x & 0x40000000))
862 ULong *
x, *xa, *xae, *xb, *xbe, *xc, *xc0;
873 if (a->wds < b->wds) {
885 for(
x = c->
x, xa =
x + wc;
x < xa;
x++)
893 for(; xb < xbe; xc0++) {
899 z = *
x++ * (ULLong)y + *xc + carry;
909 for(; xb < xbe; xb++, xc0++) {
910 if (y = *xb & 0xffff) {
915 z = (*
x & 0xffff) * y + (*xc & 0xffff) + carry;
917 z2 = (*
x++ >> 16) * y + (*xc >> 16) + carry;
930 z = (*
x & 0xffff) * y + (*xc >> 16) + carry;
933 z2 = (*
x++ >> 16) * y + (*xc & 0xffff) + carry;
941 for(; xb < xbe; xc0++) {
947 z = *
x++ * y + *xc + carry;
957 for(xc0 = c->
x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
974 static int p05[3] = { 5, 25, 125 };
983 #ifdef MULTIPLE_THREADS
1003 if (!(p51 = p5->
next)) {
1004 #ifdef MULTIPLE_THREADS
1006 if (!(p51 = p5->
next)) {
1039 n1 = n + b->wds + 1;
1040 for(i = b->maxwds; n1 > i; i <<= 1)
1044 for(i = 0; i < n; i++)
1053 *x1++ = *
x <<
k | z;
1065 *x1++ = *
x <<
k & 0xffff | z;
1089 ULong *xa, *xa0, *xb, *xb0;
1095 if (i > 1 && !a->x[i-1])
1096 Bug(
"cmp called with a->x[a->wds-1] == 0");
1097 if (j > 1 && !b->x[j-1])
1098 Bug(
"cmp called with b->x[b->wds-1] == 0");
1108 return *xa < *xb ? -1 : 1;
1125 ULong *xa, *xae, *xb, *xbe, *xc;
1162 y = (ULLong)*xa++ - *xb++ - borrow;
1163 borrow = y >> 32 & (
ULong)1;
1169 borrow = y >> 32 & (
ULong)1;
1175 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1176 borrow = (y & 0x10000) >> 16;
1177 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1178 borrow = (z & 0x10000) >> 16;
1183 y = (*xa & 0xffff) - borrow;
1184 borrow = (y & 0x10000) >> 16;
1185 z = (*xa++ >> 16) - borrow;
1186 borrow = (z & 0x10000) >> 16;
1191 y = *xa++ - *xb++ - borrow;
1192 borrow = (y & 0x10000) >> 16;
1198 borrow = (y & 0x10000) >> 16;
1221 #ifndef Avoid_Underflow
1222 #ifndef Sudden_Underflow
1231 #ifndef Avoid_Underflow
1232 #ifndef Sudden_Underflow
1237 word0(&u) = 0x80000 >> L;
1243 word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1254 (a, e)
Bigint *a;
int *e;
1259 ULong *xa, *xa0, w, y, z;
1265 #define d0 word0(&d)
1266 #define d1 word1(&d)
1273 if (!y) Bug(
"zero y in b2d");
1280 w = xa > xa0 ? *--xa : 0;
1284 z = xa > xa0 ? *--xa : 0;
1286 d0 =
Exp_1 | y << k | z >> (32 -
k);
1287 y = xa > xa0 ? *--xa : 0;
1288 d1 = z << k | y >> (32 -
k);
1296 z = xa > xa0 ? *--xa : 0;
1298 w = xa > xa0 ? *--xa : 0;
1299 y = xa > xa0 ? *--xa : 0;
1303 z = xa > xa0 ? *--xa : 0;
1304 w = xa > xa0 ? *--xa : 0;
1306 d0 =
Exp_1 | y << k + 16 | z << k | w >> 16 -
k;
1307 y = xa > xa0 ? *--xa : 0;
1308 d1 = w <<
k + 16 | y <<
k;
1324 (d, e, bits)
U *d;
int *e, *bits;
1326 (
U *d,
int *e,
int *bits)
1332 #ifndef Sudden_Underflow
1353 #ifdef Sudden_Underflow
1365 x[0] = y | z << (32 -
k);
1370 #ifndef Sudden_Underflow
1373 b->
wds = (
x[1] = z) ? 2 : 1;
1378 #ifndef Sudden_Underflow
1388 x[0] = y | z << 32 -
k & 0xffff;
1389 x[1] = z >>
k - 16 & 0xffff;
1395 x[1] = y >> 16 | z << 16 -
k & 0xffff;
1396 x[2] = z >>
k & 0xffff;
1411 Bug(
"Zero passed to d2b");
1429 #ifndef Sudden_Underflow
1433 *e = (de -
Bias - (
P-1) << 2) +
k;
1436 *e = de -
Bias - (
P-1) +
k;
1439 #ifndef Sudden_Underflow
1442 *e = de -
Bias - (
P-1) + 1 +
k;
1469 k = ka - kb + 32*(a->wds - b->wds);
1471 k = ka - kb + 16*(a->wds - b->wds);
1477 dval(&da) *= 1 <<
k;
1483 dval(&db) *= 1 <<
k;
1498 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1499 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1510 #ifdef Avoid_Underflow
1511 9007199254740992.*9007199254740992.e-256
1519 #define Scale_Bit 0x10
1523 bigtens[] = { 1e16, 1e32, 1e64 };
1548 static unsigned char hexdig[256];
1551 htinit(
unsigned char *h,
unsigned char *s,
int inc)
1554 for(i = 0; (j = s[i]) !=0; i++)
1562 #define USC (unsigned char *)
1563 htinit(
hexdig, USC
"0123456789", 0x10);
1564 htinit(
hexdig, USC
"abcdef", 0x10 + 10);
1565 htinit(
hexdig, USC
"ABCDEF", 0x10 + 10);
1569 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1570 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1571 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1572 16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
1573 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1574 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1575 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1576 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1577 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1578 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1579 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1580 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1581 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1582 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1583 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1584 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
1592 #define NAN_WORD0 0x7ff80000
1602 (sp, t)
char **sp, *t;
1604 (
const char **sp,
const char *t)
1608 CONST char *s = *sp;
1611 if ((c = *++s) >=
'A' && c <=
'Z')
1624 (rvp, sp)
U *rvp;
CONST char **sp;
1626 (
U *rvp,
const char **sp)
1631 int c1, havedig, udx0, xshift;
1635 havedig = xshift = 0;
1639 while((c = *(
CONST unsigned char*)(s+1)) && c <=
' ')
1641 if (s[1] ==
'0' && (s[2] ==
'x' || s[2] ==
'X'))
1643 while((c = *(
CONST unsigned char*)++s)) {
1646 else if (c <=
' ') {
1647 if (udx0 && havedig) {
1653 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1654 else if ( c ==
')' && havedig) {
1667 }
while((c = *++s));
1678 x[0] = (
x[0] << 4) | (
x[1] >> 28);
1679 x[1] = (
x[1] << 4) | c;
1681 if ((
x[0] &= 0xfffff) ||
x[1]) {
1699 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS)
1713 if (*
x < (
ULong)0xffffffffL) {
1754 *x1++ = (y | (*
x << n)) & 0xffffffff;
1764 if ((b->
wds = x1 - b->
x) == 0)
1783 else if (n < nwds && (
k &=
kmask)) {
1814 CONST unsigned char *decpt, *s0, *s, *s1;
1817 int big, denorm, esign, havedig,
k, n, nbits, up, zret;
1823 emax = 0x7fe -
Bias -
P + 1,
1828 emax = 0x7ff -
Bias -
P + 1
1831 emax = 0x7f -
Bias -
P
1837 #ifdef NO_LOCALE_CACHE
1838 const unsigned char *decimalpoint = (
unsigned char*)
1839 localeconv()->decimal_point;
1841 const unsigned char *decimalpoint;
1842 static unsigned char *decimalpoint_cache;
1843 if (!(s0 = decimalpoint_cache)) {
1844 s0 = (
unsigned char*)localeconv()->decimal_point;
1845 if ((decimalpoint_cache = (
unsigned char*)
1847 strcpy((
char*)decimalpoint_cache, (
CONST char*)s0);
1848 s0 = decimalpoint_cache;
1857 s0 = *(
CONST unsigned char **)sp + 2;
1858 while(s0[havedig] ==
'0')
1870 for(i = 0; decimalpoint[i]; ++i) {
1871 if (s[i] != decimalpoint[i])
1892 if (*s == *decimalpoint && !decpt) {
1893 for(i = 1; decimalpoint[i]; ++i) {
1894 if (s[i] != decimalpoint[i])
1899 if (*s ==
'.' && !decpt) {
1906 e = -(((
Long)(s-decpt)) << 2);
1920 if ((n =
hexdig[*s]) == 0 || n > 0x19) {
1925 while((n =
hexdig[*++s]) !=0 && n <= 0x19) {
1926 if (e1 & 0xf8000000)
1928 e1 = 10*e1 + n - 0x10;
1936 *sp = (
char*)s0 - 1;
1984 for(
k = 0; n > (1 << (
kshift-2)) - 1; n >>= 1)
1991 for(i = 0; decimalpoint[i+1]; ++i);
1995 if (*--s1 == decimalpoint[i]) {
2008 L |= (
hexdig[*s1] & 0x0f) << n;
2012 b->
wds = n =
x - b->
x;
2031 else if (n < nbits) {
2056 if (n == nbits && (n < 2 ||
any_on(b,n-1)))
2095 && (lostbits & 1) | (
x[0] & 1))
2110 if (nbits ==
Nbits - 1
2116 || ((n = nbits &
kmask) !=0
2126 word0(rvp) = b->
wds > 1 ? b->
x[1] & ~0x100000 : 0;
2128 word0(rvp) = (b->
x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2133 k = b->
x[0] & ((1 << j) - 1);
2147 if (
k & j && ((
k & (j-1)) | lostbits))
2153 word0(rvp) = b->
x[1] | ((e + 65 + 13) << 24);
2160 word0(rvp) = ((b->
x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->
x[1] << 16);
2161 word1(rvp) = (b->
x[0] >> 16) | (b->
x[0] << 16);
2189 ULong *bx, *bxe, q, *sx, *sxe;
2191 ULLong borrow, carry, y, ys;
2193 ULong borrow, carry, y, ys;
2202 Bug(
"oversize b in quorem");
2210 q = *bxe / (*sxe + 1);
2212 #ifdef NO_STRTOD_BIGCOMP
2219 Bug(
"oversized quotient in quorem");
2226 ys = *sx++ * (ULLong)q + carry;
2228 y = *bx - (ys &
FFFFFFFF) - borrow;
2229 borrow = y >> 32 & (
ULong)1;
2234 ys = (si & 0xffff) * q + carry;
2235 zs = (si >> 16) * q + (ys >> 16);
2237 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2238 borrow = (y & 0x10000) >> 16;
2239 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2240 borrow = (z & 0x10000) >> 16;
2243 ys = *sx++ * q + carry;
2245 y = *bx - (ys & 0xffff) - borrow;
2246 borrow = (y & 0x10000) >> 16;
2254 while(--bxe > bx && !*bxe)
2259 if (
cmp(b, S) >= 0) {
2269 y = *bx - (ys &
FFFFFFFF) - borrow;
2270 borrow = y >> 32 & (
ULong)1;
2275 ys = (si & 0xffff) + carry;
2276 zs = (si >> 16) + (ys >> 16);
2278 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2279 borrow = (y & 0x10000) >> 16;
2280 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2281 borrow = (z & 0x10000) >> 16;
2286 y = *bx - (ys & 0xffff) - borrow;
2287 borrow = (y & 0x10000) >> 16;
2296 while(--bxe > bx && !*bxe)
2304 #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP)
2326 #ifndef NO_STRTOD_BIGCOMP
2333 (
U *rv,
const char *s0,
BCinfo *bc)
2337 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2342 p5 = nd + bc->e0 - 1;
2344 #ifndef Sudden_Underflow
2350 #ifdef Avoid_Underflow
2356 #ifdef Honor_FLT_ROUNDS
2357 if (bc->rounding == 1)
2368 b =
d2b(rv, &p2, &bbits);
2369 #ifdef Avoid_Underflow
2375 if (i > (j =
P -
Emin - 1 + p2)) {
2376 #ifdef Sudden_Underflow
2381 #ifdef Avoid_Underflow
2391 #ifdef Honor_FLT_ROUNDS
2392 if (bc->rounding != 1) {
2404 #ifndef Sudden_Underflow
2433 if (!(dig =
quorem(b,d))) {
2440 for(i = 0; i < nd0; ) {
2441 if ((dd = s0[i++] -
'0' - dig))
2443 if (!b->
x[0] && b->
wds == 1) {
2451 for(j = bc->dp1; i++ < nd;) {
2452 if ((dd = s0[j++] -
'0' - dig))
2454 if (!b->
x[0] && b->
wds == 1) {
2462 if (dig > 0 || b->
x[0] || b->
wds > 1)
2467 #ifdef Honor_FLT_ROUNDS
2468 if (bc->rounding != 1) {
2470 if (bc->rounding == 0) {
2478 if (bc->rounding == 0) {
2515 if (
word1(rv) & (0x1 << i))
2518 else if (
word0(rv) & (0x1 << (i-32)))
2521 else if (
word1(rv) & 1) {
2529 #ifdef Honor_FLT_ROUNDS
2539 (s00, se)
CONST char *s00;
char **se;
2541 (
const char *s00,
char **se)
2544 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2545 int esign, i, j,
k, nd, nd0, nf, nz, nz0, nz1,
sign;
2546 CONST char *s, *s0, *s1;
2549 U aadj2, adj, rv, rv0;
2552 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2553 #ifdef Avoid_Underflow
2559 #ifndef NO_STRTOD_BIGCOMP
2560 int req_bigcomp = 0;
2562 #ifdef Honor_FLT_ROUNDS
2563 #ifdef Trust_FLT_ROUNDS
2567 switch(fegetround()) {
2568 case FE_TOWARDZERO: bc.
rounding = 0;
break;
2569 case FE_UPWARD: bc.
rounding = 2;
break;
2580 for(s = s00;;s++)
switch(*s) {
2606 #ifdef Honor_FLT_ROUNDS
2615 while(*++s ==
'0') ;
2621 for(nd = nf = 0; (c = *s) >=
'0' && c <=
'9'; nd++, s++)
2627 bc.
dp0 = bc.
dp1 = s - s0;
2628 for(s1 = s; s1 > s0 && *--s1 ==
'0'; )
2631 s1 = localeconv()->decimal_point;
2654 for(; c ==
'0'; c = *++s)
2656 if (c >
'0' && c <=
'9') {
2666 for(; c >=
'0' && c <=
'9'; c = *++s) {
2671 for(i = 1; i < nz; i++)
2674 else if (nd <= DBL_DIG + 1)
2678 else if (nd <= DBL_DIG + 1)
2686 if (c ==
'e' || c ==
'E') {
2687 if (!nd && !nz && !nz0) {
2698 if (c >=
'0' && c <=
'9') {
2701 if (c >
'0' && c <=
'9') {
2704 while((c = *++s) >=
'0' && c <=
'9')
2706 if (s - s1 > 8 || L > 19999)
2730 if (
match(&s,
"nf")) {
2732 if (!
match(&s,
"inity"))
2734 word0(&rv) = 0x7ff00000;
2741 if (
match(&s,
"an")) {
2758 bc.
e0 = e1 = e -= nf;
2767 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2772 oldinexact = get_inexact();
2778 #ifndef RND_PRODQUOT
2779 #ifndef Honor_FLT_ROUNDS
2786 #ifndef ROUND_BIASED_without_Round_Up
2790 goto vax_ovfl_check;
2792 #ifdef Honor_FLT_ROUNDS
2808 #ifdef Honor_FLT_ROUNDS
2834 #ifndef Inaccurate_Divide
2836 #ifdef Honor_FLT_ROUNDS
2855 oldinexact = get_inexact();
2857 #ifdef Avoid_Underflow
2860 #ifdef Honor_FLT_ROUNDS
2877 if (e1 > DBL_MAX_10_EXP) {
2881 #ifdef Honor_FLT_ROUNDS
2919 for(j = 0; e1 > 1; j++, e1 >>= 1)
2945 #ifdef Avoid_Underflow
2948 for(j = 0; e1 > 0; j++, e1 >>= 1)
2961 word0(&rv) &= 0xffffffff << (j-32);
2964 word1(&rv) &= 0xffffffff << j;
2967 for(j = 0; e1 > 1; j++, e1 >>= 1)
2982 #ifndef Avoid_Underflow
2998 #ifndef NO_STRTOD_BIGCOMP
3010 if (--j < bc.dp1 && j >= bc.
dp0)
3022 for(i = 0; i < nd0; ++i)
3023 y = 10*y + s0[i] -
'0';
3024 for(j = bc.
dp1; i < nd; ++i)
3025 y = 10*y + s0[j++] -
'0';
3029 bd0 =
s2b(s0, nd0, nd, y, bc.
dplen);
3034 bb =
d2b(&rv, &bbe, &bbbits);
3050 #ifdef Honor_FLT_ROUNDS
3054 #ifdef Avoid_Underflow
3066 Lsb1 = Lsb << (i-32);
3071 #ifdef Sudden_Underflow
3073 j = 1 + 4*
P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
3088 #ifdef Avoid_Underflow
3091 i = bb2 < bd2 ? bb2 : bd2;
3113 delta =
diff(bb, bd);
3117 #ifndef NO_STRTOD_BIGCOMP
3118 if (bc.
nd > nd && i <= 0) {
3124 #ifdef Honor_FLT_ROUNDS
3136 #ifdef Honor_FLT_ROUNDS
3140 if (!delta->
x[0] && delta->
wds <= 1) {
3153 else if (!bc.
dsign) {
3158 #ifdef Avoid_Underflow
3165 if (
cmp(delta, bs) <= 0)
3170 #ifdef Avoid_Underflow
3175 #ifdef Sudden_Underflow
3189 adj.
d =
ratio(delta, bs);
3192 if (adj.
d <= 0x7ffffffe) {
3201 #ifdef Avoid_Underflow
3205 #ifdef Sudden_Underflow
3244 if (!delta->
x[0] && delta->
wds <= 1)
3249 if (!delta->
x[0] && delta->
wds <= 1) {
3257 if (
cmp(delta, bs) > 0)
3268 ? (0xffffffff & (0xffffffff << (2*
P+1-(y>>
Exp_shift)))) :
3281 #ifdef Avoid_Underflow
3290 #ifdef Sudden_Underflow
3295 #ifdef Avoid_Underflow
3310 #ifdef Avoid_Underflow
3330 word1(&rv) = 0xffffffff;
3334 #ifndef NO_STRTOD_BIGCOMP
3341 #ifndef ROUND_BIASED
3342 #ifdef Avoid_Underflow
3344 if (!(
word0(&rv) & Lsb1))
3347 else if (!(
word1(&rv) & Lsb))
3355 #ifdef Avoid_Underflow
3360 #ifndef ROUND_BIASED
3362 #ifdef Avoid_Underflow
3367 #ifndef Sudden_Underflow
3377 #ifdef Avoid_Underflow
3383 if ((aadj =
ratio(delta, bs)) <= 2.) {
3387 #ifndef Sudden_Underflow
3403 if (aadj < 2./FLT_RADIX)
3404 aadj = 1./FLT_RADIX;
3412 aadj1 = bc.
dsign ? aadj : -aadj;
3413 #ifdef Check_FLT_ROUNDS
3434 adj.
d = aadj1 *
ulp(&rv);
3448 #ifdef Avoid_Underflow
3450 if (aadj <= 0x7fffffff) {
3451 if ((z = aadj) <= 0)
3454 aadj1 = bc.
dsign ? aadj : -aadj;
3456 dval(&aadj2) = aadj1;
3458 aadj1 =
dval(&aadj2);
3459 adj.
d = aadj1 *
ulp(&rv);
3462 #ifdef NO_STRTOD_BIGCOMP
3473 adj.
d = aadj1 *
ulp(&rv);
3477 #ifdef Sudden_Underflow
3481 adj.
d = aadj1 *
ulp(&rv);
3505 adj.
d = aadj1 *
ulp(&rv);
3517 aadj1 = (double)(
int)(aadj + 0.5);
3521 adj.
d = aadj1 *
ulp(&rv);
3529 #ifdef Avoid_Underflow
3538 if (aadj < .4999999 || aadj > .5000001)
3541 else if (aadj < .4999999/FLT_RADIX)
3557 #ifndef NO_STRTOD_BIGCOMP
3565 if (y == 0 && rv.
d == 0.)
3577 else if (!oldinexact)
3580 #ifdef Avoid_Underflow
3599 dval(&rv0) = 1e-300;
3609 #ifndef MULTIPLE_THREADS
3624 (int)(
sizeof(
Bigint) -
sizeof(
ULong) -
sizeof(
int)) + j <= i;
3630 #ifndef MULTIPLE_THREADS
3638 nrv_alloc(s, rve, n)
char *s, **rve;
int n;
3646 while((*t = *s++)) t++;
3666 b->
maxwds = 1 << (b->
k = *(
int*)b);
3668 #ifndef MULTIPLE_THREADS
3711 (dd, mode, ndigits, decpt,
sign, rve)
3712 double dd;
int mode, ndigits, *decpt, *
sign;
char **rve;
3714 (
double dd,
int mode,
int ndigits,
int *decpt,
int *
sign,
char **rve)
3751 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3752 j, j1,
k, k0, k_check, leftright, m2, m5, s2, s5,
3753 spec_case, try_quick;
3755 #ifndef Sudden_Underflow
3759 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3763 #ifndef No_leftright
3769 int inexact, oldinexact;
3771 #ifdef Honor_FLT_ROUNDS
3773 #ifdef Trust_FLT_ROUNDS
3777 switch(fegetround()) {
3778 case FE_TOWARDZERO:
Rounding = 0;
break;
3779 case FE_UPWARD:
Rounding = 2;
break;
3785 #ifndef MULTIPLE_THREADS
3801 #if defined(IEEE_Arith) + defined(VAX)
3805 if (
word0(&u) == 0x8000)
3826 try_quick = oldinexact = get_inexact();
3829 #ifdef Honor_FLT_ROUNDS
3839 b =
d2b(&u, &be, &bbits);
3840 #ifdef Sudden_Underflow
3850 dval(&d2) /= 1 << j;
3880 #ifndef Sudden_Underflow
3886 i = bbits + be + (
Bias + (
P-1) - 1);
3887 x = i > 32 ?
word0(&u) << (64 - i) |
word1(&u) >> (i - 32)
3888 :
word1(&u) << (32 - i);
3891 i -= (
Bias + (
P-1) - 1) + 1;
3895 ds = (
dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3897 if (ds < 0. && ds !=
k)
3924 if (mode < 0 || mode > 9)
3928 #ifdef Check_FLT_ROUNDS
3954 ilim = ilim1 = i = ndigits;
3960 i = ndigits +
k + 1;
3968 #ifdef Honor_FLT_ROUNDS
3973 if (ilim >= 0 && ilim <=
Quick_max && try_quick) {
3992 for(; j; j >>= 1, i++)
3999 else if ((j1 = -
k)) {
4001 for(j = j1 >> 4; j; j >>= 1, i++)
4007 if (k_check &&
dval(&u) < 1. && ilim > 0) {
4026 #ifndef No_leftright
4033 if (k0 < 0 && j1 >= 307) {
4037 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
4047 *s++ =
'0' + (int)L;
4062 for(i = 1;; i++,
dval(&u) *= 10.) {
4064 if (!(
dval(&u) -= L))
4066 *s++ =
'0' + (int)L;
4070 else if (
dval(&u) < 0.5 -
dval(&eps)) {
4078 #ifndef No_leftright
4093 if (ndigits < 0 && ilim <= 0) {
4095 if (ilim < 0 ||
dval(&u) <= 5*ds)
4099 for(i = 1;; i++,
dval(&u) *= 10.) {
4102 #ifdef Check_FLT_ROUNDS
4109 *s++ =
'0' + (int)L;
4117 #ifdef Honor_FLT_ROUNDS
4121 case 2:
goto bump_up;
4128 if (
dval(&u) > ds || (
dval(&u) == ds && L & 1))
4151 #ifndef Sudden_Underflow
4152 denorm ? be + (
Bias + (
P-1) - 1 + 1) :
4155 1 + 4*
P - 3 - bbits + ((bbits + be - 1) & 3);
4163 if (m2 > 0 && s2 > 0) {
4164 i = m2 < s2 ? m2 : s2;
4190 if ((mode < 2 || leftright)
4191 #ifdef Honor_FLT_ROUNDS
4196 #ifndef Sudden_Underflow
4231 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4232 if (ilim < 0 ||
cmp(b,S =
multadd(S,5,0)) <= 0) {
4264 delta =
diff(S, mhi);
4265 j1 = delta->
sign ? 1 :
cmp(b, delta);
4267 #ifndef ROUND_BIASED
4268 if (j1 == 0 && mode != 1 && !(
word1(&u) & 1)
4269 #ifdef Honor_FLT_ROUNDS
4278 else if (!b->
x[0] && b->
wds <= 1)
4285 if (j < 0 || (j == 0 && mode != 1
4286 #ifndef ROUND_BIASED
4290 if (!b->
x[0] && b->
wds <= 1) {
4296 #ifdef Honor_FLT_ROUNDS
4299 case 0:
goto accept_dig;
4300 case 2:
goto keep_dig;
4309 if ((j1 > 0 || (j1 == 0 && dig & 1))
4319 #ifdef Honor_FLT_ROUNDS
4331 #ifdef Honor_FLT_ROUNDS
4339 mlo = mhi =
multadd(mhi, 10, 0);
4348 *s++ = dig =
quorem(b,S) +
'0';
4349 if (!b->
x[0] && b->
wds <= 1) {
4362 #ifdef Honor_FLT_ROUNDS
4364 case 0:
goto trimzeros;
4365 case 2:
goto roundoff;
4373 if (j > 0 || (j == 0 && dig & 1))
4386 #ifdef Honor_FLT_ROUNDS
4395 if (mlo && mlo != mhi)
4408 else if (!oldinexact)
#define rounded_product(a, b)
static Bigint * pow5mult(Bigint *b, int k)
static void bigcomp(U *rv, const char *s0, BCinfo *bc)
static char * rv_alloc(int i)
static double b2d(Bigint *a, int *e)
static Bigint * d2b(U *d, int *e, int *bits)
#define FREE_DTOA_LOCK(n)
static int lo0bits(ULong *y)
static void rshift(Bigint *b, int k)
#define rounded_quotient(a, b)
static Bigint * increment(Bigint *b)
double os_strtod(const char *s00, char **se)
static int quorem(Bigint *b, Bigint *S)
void os_freedtoa(char *s)
static CONST double tinytens[]
static int match(const char **sp, const char *t)
static Bigint * s2b(const char *s, int nd0, int nd, ULong y9, int dplen)
static Bigint * i2b(int i)
static Bigint * diff(Bigint *a, Bigint *b)
static void Bfree(Bigint *v)
static Bigint * multadd(Bigint *b, int m, int a)
void os_gethex(CONST char **sp, U *rvp, int rounding, int sign)
static Bigint * lshift(Bigint *b, int k)
#define Storeinc(a, b, c)
static Bigint * Balloc(int k)
#define ACQUIRE_DTOA_LOCK(n)
static int dshift(Bigint *b, int p2)
static void hexnan(U *rvp, const char **sp)
char * os_dtoa(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
static Bigint * freelist[Kmax+1]
static char * dtoa_result
static int hi0bits(ULong x)
static ULong any_on(Bigint *b, int k)
static double sulp(U *x, BCinfo *bc)
static int cmp(Bigint *a, Bigint *b)
static double ratio(Bigint *a, Bigint *b)
static CONST double tens[]
static double private_mem[PRIVATE_mem]
static Bigint * mult(Bigint *a, Bigint *b)
static CONST double bigtens[]
static char * nrv_alloc(const char *s, char **rve, int n)
static double * pmem_next
unsigned Long ULong
end of OS code, below is David Gay, except as noted above
static unsigned char hexdig[256]