166 #ifdef WORDS_BIGENDIAN
167 #define IEEE_BIG_ENDIAN
169 #define IEEE_LITTLE_ENDIAN
174 #undef IEEE_BIG_ENDIAN
175 #undef IEEE_LITTLE_ENDIAN
178 #if defined(__arm__) && !defined(__VFP_FP__)
179 #define IEEE_BIG_ENDIAN
180 #undef IEEE_LITTLE_ENDIAN
191 #if (INT_MAX >> 30) && !(INT_MAX >> 31)
193 #define ULong unsigned int
194 #elif (LONG_MAX >> 30) && !(LONG_MAX >> 31)
195 #define Long long int
196 #define ULong unsigned long int
198 #error No 32bit integer
201 #if defined(HAVE_LONG_LONG) && (HAVE_LONG_LONG)
202 #define Llong LONG_LONG
209 #define Bug(x) {fprintf(stderr, "%s\n", (x)); exit(EXIT_FAILURE);}
214 #define ISDIGIT(c) isdigit(c)
225 extern void *MALLOC(
size_t);
227 #define MALLOC malloc
230 extern void FREE(
void*);
235 #define NO_SANITIZE(x, y) y
238 #ifndef Omit_Private_Memory
240 #define PRIVATE_MEM 2304
242 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
243 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
247 #undef Avoid_Underflow
248 #ifdef IEEE_BIG_ENDIAN
251 #ifdef IEEE_LITTLE_ENDIAN
259 #define DBL_MAX_10_EXP 308
260 #define DBL_MAX_EXP 1024
266 #define DBL_MAX_10_EXP 75
267 #define DBL_MAX_EXP 63
269 #define DBL_MAX 7.2370055773322621e+75
274 #define DBL_MAX_10_EXP 38
275 #define DBL_MAX_EXP 127
277 #define DBL_MAX 1.7014118346046923e+38
281 #define LONG_MAX 2147483647
298 static const char hexdigit[] =
"0123456789abcdef0123456789ABCDEF";
301 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + defined(IBM) != 1
302 Exactly one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM should be defined.
305 typedef union {
double d; ULong L[2]; }
U;
310 # ifdef IEEE_LITTLE_ENDIAN
311 # define word0(x) (((ULong *)&(x))[1])
312 # define word1(x) (((ULong *)&(x))[0])
314 # define word0(x) (((ULong *)&(x))[0])
315 # define word1(x) (((ULong *)&(x))[1])
319 # ifdef IEEE_LITTLE_ENDIAN
320 # define word0(x) ((x).L[1])
321 # define word1(x) ((x).L[0])
323 # define word0(x) ((x).L[0])
324 # define word1(x) ((x).L[1])
326 # define dval(x) ((x).d)
333 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
334 #define Storeinc(a,b,c) (((unsigned short *)(a))[1] = (unsigned short)(b), \
335 ((unsigned short *)(a))[0] = (unsigned short)(c), (a)++)
337 #define Storeinc(a,b,c) (((unsigned short *)(a))[0] = (unsigned short)(b), \
338 ((unsigned short *)(a))[1] = (unsigned short)(c), (a)++)
349 #define Exp_shift1 20
350 #define Exp_msk1 0x100000
351 #define Exp_msk11 0x100000
352 #define Exp_mask 0x7ff00000
356 #define Exp_1 0x3ff00000
357 #define Exp_11 0x3ff00000
359 #define Frac_mask 0xfffff
360 #define Frac_mask1 0xfffff
363 #define Bndry_mask 0xfffff
364 #define Bndry_mask1 0xfffff
366 #define Sign_bit 0x80000000
372 #ifndef NO_IEEE_Scale
373 #define Avoid_Underflow
375 #undef Sudden_Underflow
381 #define Flt_Rounds FLT_ROUNDS
387 #ifdef Honor_FLT_ROUNDS
388 #define Rounding rounding
389 #undef Check_FLT_ROUNDS
390 #define Check_FLT_ROUNDS
392 #define Rounding Flt_Rounds
396 #undef Check_FLT_ROUNDS
397 #undef Honor_FLT_ROUNDS
399 #undef Sudden_Underflow
400 #define Sudden_Underflow
405 #define Exp_shift1 24
406 #define Exp_msk1 0x1000000
407 #define Exp_msk11 0x1000000
408 #define Exp_mask 0x7f000000
411 #define Exp_1 0x41000000
412 #define Exp_11 0x41000000
414 #define Frac_mask 0xffffff
415 #define Frac_mask1 0xffffff
418 #define Bndry_mask 0xefffff
419 #define Bndry_mask1 0xffffff
421 #define Sign_bit 0x80000000
423 #define Tiny0 0x100000
432 #define Exp_msk1 0x80
433 #define Exp_msk11 0x800000
434 #define Exp_mask 0x7f80
437 #define Exp_1 0x40800000
438 #define Exp_11 0x4080
440 #define Frac_mask 0x7fffff
441 #define Frac_mask1 0xffff007f
444 #define Bndry_mask 0xffff007f
445 #define Bndry_mask1 0xffff007f
447 #define Sign_bit 0x8000
461 #define rounded_product(a,b) ((a) = rnd_prod((a), (b)))
462 #define rounded_quotient(a,b) ((a) = rnd_quot((a), (b)))
463 extern double rnd_prod(
double,
double), rnd_quot(
double,
double);
465 #define rounded_product(a,b) ((a) *= (b))
466 #define rounded_quotient(a,b) ((a) /= (b))
469 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
470 #define Big1 0xffffffff
476 #define FFFFFFFF 0xffffffffUL
490 #define Llong long long
493 #define ULLong unsigned Llong
497 #define MULTIPLE_THREADS 1
499 #ifndef MULTIPLE_THREADS
500 #define ACQUIRE_DTOA_LOCK(n)
501 #define FREE_DTOA_LOCK(n)
503 #define ACQUIRE_DTOA_LOCK(n)
504 #define FREE_DTOA_LOCK(n)
507 #ifndef ATOMIC_PTR_CAS
508 #define ATOMIC_PTR_CAS(var, old, new) ((var) = (new), (void *)(old))
511 #define LIKELY(x) (x)
514 #define UNLIKELY(x) (x)
517 #define ASSUME(x) (void)(x)
524 int k, maxwds, sign, wds;
530 static Bigint *freelist[Kmax+1];
532 #define BLOCKING_BIGINT ((Bigint *)(-1))
539 #ifndef Omit_Private_Memory
544 ACQUIRE_DTOA_LOCK(0);
549 rv = ATOMIC_PTR_CAS(freelist[k], rv, BLOCKING_BIGINT);
550 if (LIKELY(rv != BLOCKING_BIGINT && rvn == rv)) {
551 rvn = ATOMIC_PTR_CAS(freelist[k], BLOCKING_BIGINT, rv->next);
552 assert(rvn == BLOCKING_BIGINT);
560 #ifdef Omit_Private_Memory
561 rv = (
Bigint *)MALLOC(
sizeof(
Bigint) + (x-1)*
sizeof(ULong));
563 len = (
sizeof(
Bigint) + (x-1)*
sizeof(ULong) +
sizeof(
double) - 1)
566 double *pnext = pmem_next;
567 while (pnext - private_mem +
len <= PRIVATE_mem) {
569 pnext = ATOMIC_PTR_CAS(pmem_next, pnext, pnext +
len);
570 if (LIKELY(p == pnext)) {
578 rv = (
Bigint*)MALLOC(
len*
sizeof(
double));
584 rv->sign = rv->wds = 0;
597 ACQUIRE_DTOA_LOCK(0);
600 vn = ATOMIC_PTR_CAS(freelist[v->k], 0, 0);
601 }
while (UNLIKELY(vn == BLOCKING_BIGINT));
603 }
while (UNLIKELY(ATOMIC_PTR_CAS(freelist[v->k], vn, v) != vn));
608 #define Bcopy(x,y) memcpy((char *)&(x)->sign, (char *)&(y)->sign, \
609 (y)->wds*sizeof(Long) + 2*sizeof(int))
612 multadd(
Bigint *b,
int m,
int a)
632 y = *x * (ULLong)m + carry;
634 *x++ = (ULong)(y & FFFFFFFF);
638 y = (xi & 0xffff) * m + carry;
639 z = (xi >> 16) * m + (y >> 16);
641 *x++ = (z << 16) + (y & 0xffff);
650 if (wds >= b->maxwds) {
656 b->x[wds++] = (ULong)carry;
663 s2b(
const char *s,
int nd0,
int nd, ULong y9)
670 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
677 b->x[0] = y9 & 0xffff;
678 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
685 b = multadd(b, 10, *s++ -
'0');
692 b = multadd(b, 10, *s++ -
'0');
697 hi0bits(
register ULong x)
701 if (!(x & 0xffff0000)) {
705 if (!(x & 0xff000000)) {
709 if (!(x & 0xf0000000)) {
713 if (!(x & 0xc0000000)) {
717 if (!(x & 0x80000000)) {
719 if (!(x & 0x40000000))
729 register ULong x = *y;
784 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
795 if (a->wds < b->wds) {
807 for (x = c->x, xa = x + wc; x < xa; x++)
815 for (; xb < xbe; xc0++) {
816 if ((y = *xb++) != 0) {
821 z = *x++ * (ULLong)y + *xc + carry;
823 *xc++ = (ULong)(z & FFFFFFFF);
830 for (; xb < xbe; xb++, xc0++) {
831 if ((y = *xb & 0xffff) != 0) {
836 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
838 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
844 if ((y = *xb >> 16) != 0) {
850 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
853 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
860 for (; xb < xbe; xc0++) {
866 z = *x++ * y + *xc + carry;
875 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
883 pow5mult(
Bigint *b,
int k)
888 static const int p05[3] = { 5, 25, 125 };
890 if ((i = k & 3) != 0)
891 b = multadd(b, p05[i-1], 0);
897 ACQUIRE_DTOA_LOCK(1);
901 p5tmp = ATOMIC_PTR_CAS(p5s, NULL, p5);
902 if (UNLIKELY(p5tmp)) {
917 if (!(p51 = p5->next)) {
918 ACQUIRE_DTOA_LOCK(1);
919 if (!(p51 = p5->next)) {
922 p5tmp = ATOMIC_PTR_CAS(p5->next, NULL, p51);
923 if (UNLIKELY(p5tmp)) {
940 ULong *x, *x1, *xe, z;
949 for (i = b->maxwds; n1 > i; i <<= 1)
953 for (i = 0; i < n; i++)
973 *x1++ = *x << k & 0xffff | z;
992 ULong *xa, *xa0, *xb, *xb0;
998 if (i > 1 && !a->x[i-1])
999 Bug(
"cmp called with a->x[a->wds-1] == 0");
1000 if (j > 1 && !b->x[j-1])
1001 Bug(
"cmp called with b->x[b->wds-1] == 0");
1011 return *xa < *xb ? -1 : 1;
1018 NO_SANITIZE(
"unsigned-integer-overflow",
static Bigint * diff(
Bigint *a,
Bigint *b));
1024 ULong *xa, *xae, *xb, *xbe, *xc;
1061 y = (ULLong)*xa++ - *xb++ - borrow;
1062 borrow = y >> 32 & (ULong)1;
1063 *xc++ = (ULong)(y & FFFFFFFF);
1067 borrow = y >> 32 & (ULong)1;
1068 *xc++ = (ULong)(y & FFFFFFFF);
1073 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1074 borrow = (y & 0x10000) >> 16;
1075 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1076 borrow = (z & 0x10000) >> 16;
1080 y = (*xa & 0xffff) - borrow;
1081 borrow = (y & 0x10000) >> 16;
1082 z = (*xa++ >> 16) - borrow;
1083 borrow = (z & 0x10000) >> 16;
1088 y = *xa++ - *xb++ - borrow;
1089 borrow = (y & 0x10000) >> 16;
1094 borrow = (y & 0x10000) >> 16;
1112 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1113 #ifndef Avoid_Underflow
1114 #ifndef Sudden_Underflow
1123 #ifndef Avoid_Underflow
1124 #ifndef Sudden_Underflow
1127 L = -L >> Exp_shift;
1128 if (L < Exp_shift) {
1129 word0(a) = 0x80000 >> L;
1135 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1146 ULong *xa, *xa0, w, y, z;
1160 if (!y) Bug(
"zero y in b2d");
1166 d0 = Exp_1 | y >> (Ebits - k);
1167 w = xa > xa0 ? *--xa : 0;
1168 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1171 z = xa > xa0 ? *--xa : 0;
1173 d0 = Exp_1 | y << k | z >> (32 - k);
1174 y = xa > xa0 ? *--xa : 0;
1175 d1 = z << k | y >> (32 - k);
1182 if (k < Ebits + 16) {
1183 z = xa > xa0 ? *--xa : 0;
1184 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1185 w = xa > xa0 ? *--xa : 0;
1186 y = xa > xa0 ? *--xa : 0;
1187 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1190 z = xa > xa0 ? *--xa : 0;
1191 w = xa > xa0 ? *--xa : 0;
1193 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1194 y = xa > xa0 ? *--xa : 0;
1195 d1 = w << k + 16 | y << k;
1199 word0(d) = d0 >> 16 | d0 << 16;
1200 word1(d) = d1 >> 16 | d1 << 16;
1209 d2b(
double d_,
int *e,
int *bits)
1215 #ifndef Sudden_Underflow
1223 d0 = word0(d) >> 16 | word0(d) << 16;
1224 d1 = word1(d) >> 16 | word1(d) << 16;
1239 #ifdef Sudden_Underflow
1240 de = (int)(d0 >> Exp_shift);
1245 if ((de = (
int)(d0 >> Exp_shift)) != 0)
1249 if ((y = d1) != 0) {
1250 if ((k = lo0bits(&y)) != 0) {
1251 x[0] = y | z << (32 - k);
1256 #ifndef Sudden_Underflow
1259 b->wds = (x[1] = z) ? 2 : 1;
1264 Bug(
"Zero passed to d2b");
1268 #ifndef Sudden_Underflow
1276 if (k = lo0bits(&y))
1278 x[0] = y | z << 32 - k & 0xffff;
1279 x[1] = z >> k - 16 & 0xffff;
1285 x[1] = y >> 16 | z << 16 - k & 0xffff;
1286 x[2] = z >> k & 0xffff;
1301 Bug(
"Zero passed to d2b");
1319 #ifndef Sudden_Underflow
1323 *e = (de - Bias - (P-1) << 2) + k;
1324 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1326 *e = de - Bias - (P-1) + k;
1329 #ifndef Sudden_Underflow
1332 *e = de - Bias - (P-1) + 1 + k;
1334 *bits = 32*i - hi0bits(x[i-1]);
1336 *bits = (i+2)*16 - hi0bits(x[i]);
1351 dval(da) = b2d(a, &ka);
1352 dval(db) = b2d(b, &kb);
1354 k = ka - kb + 32*(a->wds - b->wds);
1356 k = ka - kb + 16*(a->wds - b->wds);
1360 word0(da) += (k >> 2)*Exp_msk1;
1366 word0(db) += (k >> 2)*Exp_msk1;
1372 word0(da) += k*Exp_msk1;
1375 word0(db) += k*Exp_msk1;
1378 return dval(da) / dval(db);
1383 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1384 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1393 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1394 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1395 #ifdef Avoid_Underflow
1396 9007199254740992.*9007199254740992.e-256
1404 #define Scale_Bit 0x10
1408 bigtens[] = { 1e16, 1e32, 1e64 };
1409 static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1412 bigtens[] = { 1e16, 1e32 };
1413 static const double tinytens[] = { 1e-16, 1e-32 };
1425 #define NAN_WORD0 0x7ff80000
1433 match(
const char **sp,
char *t)
1436 const char *s = *sp;
1439 if ((c = *++s) >=
'A' && c <=
'Z')
1450 hexnan(
double *rvp,
const char **sp)
1454 int havedig, udx0, xshift;
1457 havedig = xshift = 0;
1460 while (c = *(
const unsigned char*)++s) {
1461 if (c >=
'0' && c <=
'9')
1463 else if (c >=
'a' && c <=
'f')
1465 else if (c >=
'A' && c <=
'F')
1467 else if (c <=
' ') {
1468 if (udx0 && havedig) {
1474 else if ( c ==
')' && havedig) {
1487 x[0] = (x[0] << 4) | (x[1] >> 28);
1488 x[1] = (x[1] << 4) | c;
1490 if ((x[0] &= 0xfffff) || x[1]) {
1491 word0(*rvp) = Exp_mask | x[0];
1498 NO_SANITIZE(
"unsigned-integer-overflow",
double strtod(
const char *s00,
char **se));
1500 strtod(
const char *s00,
char **se)
1502 #ifdef Avoid_Underflow
1505 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1506 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1507 const char *s, *s0, *s1;
1512 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1514 int inexact, oldinexact;
1516 #ifdef Honor_FLT_ROUNDS
1524 sign = nz0 = nz = 0;
1549 if (s[1] ==
'x' || s[1] ==
'X') {
1555 if (!*++s || (!(s1 = strchr(hexdigit, *s)) && *s !=
'.'))
goto ret0;
1557 while (*++s ==
'0');
1559 s1 = strchr(hexdigit, *s);
1563 adj += aadj * ((s1 - hexdigit) & 15);
1566 }
while (*++s && (s1 = strchr(hexdigit, *s)));
1569 if ((*s ==
'.') && *++s && (s1 = strchr(hexdigit, *s))) {
1576 for (; *s && (s1 = strchr(hexdigit, *s)); ++s) {
1577 adj += aadj * ((s1 - hexdigit) & 15);
1578 if ((aadj /= 16) == 0.0) {
1579 while (*++s && strchr(hexdigit, *s));
1585 if (*s ==
'P' || *s ==
'p') {
1586 dsign = 0x2C - *++s;
1587 if (abs(dsign) == 1) s++;
1592 if (c <
'0' ||
'9' < c)
goto ret0;
1599 if (nd + dsign * nd0 > 2095) {
1600 while (
'0' <= c && c <=
'9') c = *++s;
1603 }
while (
'0' <= c && c <=
'9');
1606 dval(rv) = ldexp(adj, nd0);
1610 while (*++s ==
'0') ;
1616 for (nd = nf = 0; (c = *s) >=
'0' && c <=
'9'; nd++, s++)
1619 else if (nd < DBL_DIG + 2)
1623 s1 = localeconv()->decimal_point;
1646 for (; c ==
'0'; c = *++s)
1648 if (c >
'0' && c <=
'9') {
1656 for (; c >=
'0' && c <=
'9'; c = *++s) {
1659 if (nd > DBL_DIG * 4) {
1664 for (i = 1; i < nz; i++)
1667 else if (nd <= DBL_DIG + 2)
1671 else if (nd <= DBL_DIG + 2)
1679 if (c ==
'e' || c ==
'E') {
1680 if (!nd && !nz && !nz0) {
1691 if (c >=
'0' && c <=
'9') {
1694 if (c >
'0' && c <=
'9') {
1697 while ((c = *++s) >=
'0' && c <=
'9')
1699 if (s - s1 > 8 || L > 19999)
1722 if (match(&s,
"nf")) {
1724 if (!match(&s,
"inity"))
1726 word0(rv) = 0x7ff00000;
1733 if (match(&s,
"an")) {
1734 word0(rv) = NAN_WORD0;
1735 word1(rv) = NAN_WORD1;
1759 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
1764 oldinexact = get_inexact();
1766 dval(rv) = tens[k - 9] * dval(rv) + z;
1768 bd0 = bb = bd = bs = delta = 0;
1770 #ifndef RND_PRODQUOT
1771 #ifndef Honor_FLT_ROUNDS
1779 if (e <= Ten_pmax) {
1781 goto vax_ovfl_check;
1783 #ifdef Honor_FLT_ROUNDS
1786 dval(rv) = -dval(rv);
1790 rounded_product(dval(rv), tens[e]);
1795 if (e <= Ten_pmax + i) {
1799 #ifdef Honor_FLT_ROUNDS
1802 dval(rv) = -dval(rv);
1807 dval(rv) *= tens[i];
1813 word0(rv) -= P*Exp_msk1;
1814 rounded_product(dval(rv), tens[e]);
1815 if ((word0(rv) & Exp_mask)
1816 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1818 word0(rv) += P*Exp_msk1;
1820 rounded_product(dval(rv), tens[e]);
1825 #ifndef Inaccurate_Divide
1826 else if (e >= -Ten_pmax) {
1827 #ifdef Honor_FLT_ROUNDS
1830 dval(rv) = -dval(rv);
1834 rounded_quotient(dval(rv), tens[-e]);
1845 oldinexact = get_inexact();
1847 #ifdef Avoid_Underflow
1850 #ifdef Honor_FLT_ROUNDS
1851 if ((rounding = Flt_Rounds) >= 2) {
1853 rounding = rounding == 2 ? 0 : 2;
1864 if ((i = e1 & 15) != 0)
1865 dval(rv) *= tens[i];
1867 if (e1 > DBL_MAX_10_EXP) {
1874 #ifdef Honor_FLT_ROUNDS
1882 word0(rv) = Exp_mask;
1886 word0(rv) = Exp_mask;
1892 dval(rv0) *= dval(rv0);
1903 for (j = 0; e1 > 1; j++, e1 >>= 1)
1905 dval(rv) *= bigtens[j];
1907 word0(rv) -= P*Exp_msk1;
1908 dval(rv) *= bigtens[j];
1909 if ((z = word0(rv) & Exp_mask)
1910 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1912 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1919 word0(rv) += P*Exp_msk1;
1924 if ((i = e1 & 15) != 0)
1925 dval(rv) /= tens[i];
1927 if (e1 >= 1 << n_bigtens)
1929 #ifdef Avoid_Underflow
1932 for (j = 0; e1 > 0; j++, e1 >>= 1)
1934 dval(rv) *= tinytens[j];
1935 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1936 >> Exp_shift)) > 0) {
1941 word0(rv) = (P+2)*Exp_msk1;
1943 word0(rv) &= 0xffffffff << (j-32);
1946 word1(rv) &= 0xffffffff << j;
1949 for (j = 0; e1 > 1; j++, e1 >>= 1)
1951 dval(rv) *= tinytens[j];
1953 dval(rv0) = dval(rv);
1954 dval(rv) *= tinytens[j];
1956 dval(rv) = 2.*dval(rv0);
1957 dval(rv) *= tinytens[j];
1969 #ifndef Avoid_Underflow
1984 bd0 = s2b(s0, nd0, nd, y);
1987 bd = Balloc(bd0->k);
1989 bb = d2b(dval(rv), &bbe, &bbbits);
2005 #ifdef Honor_FLT_ROUNDS
2009 #ifdef Avoid_Underflow
2017 #ifdef Sudden_Underflow
2019 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2034 #ifdef Avoid_Underflow
2037 i = bb2 < bd2 ? bb2 : bd2;
2046 bs = pow5mult(bs, bb5);
2052 bb = lshift(bb, bb2);
2054 bd = pow5mult(bd, bd5);
2056 bd = lshift(bd, bd2);
2058 bs = lshift(bs, bs2);
2059 delta = diff(bb, bd);
2060 dsign = delta->sign;
2063 #ifdef Honor_FLT_ROUNDS
2064 if (rounding != 1) {
2067 if (!delta->x[0] && delta->wds <= 1) {
2083 && !(word0(rv) & Frac_mask)) {
2084 y = word0(rv) & Exp_mask;
2085 #ifdef Avoid_Underflow
2086 if (!scale || y > 2*P*Exp_msk1)
2091 delta = lshift(delta,Log2P);
2092 if (cmp(delta, bs) <= 0)
2097 #ifdef Avoid_Underflow
2098 if (scale && (y = word0(rv) & Exp_mask)
2100 word0(adj) += (2*P+1)*Exp_msk1 - y;
2102 #ifdef Sudden_Underflow
2103 if ((word0(rv) & Exp_mask) <=
2105 word0(rv) += P*Exp_msk1;
2106 dval(rv) += adj*ulp(dval(rv));
2107 word0(rv) -= P*Exp_msk1;
2112 dval(rv) += adj*ulp(dval(rv));
2116 adj = ratio(delta, bs);
2119 if (adj <= 0x7ffffffe) {
2123 if (!((rounding>>1) ^ dsign))
2128 #ifdef Avoid_Underflow
2129 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2130 word0(adj) += (2*P+1)*Exp_msk1 - y;
2132 #ifdef Sudden_Underflow
2133 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2134 word0(rv) += P*Exp_msk1;
2135 adj *= ulp(dval(rv));
2140 word0(rv) -= P*Exp_msk1;
2145 adj *= ulp(dval(rv));
2158 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2160 #ifdef Avoid_Underflow
2161 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2163 || (word0(rv) & Exp_mask) <= Exp_msk1
2168 if (!delta->x[0] && delta->wds <= 1)
2173 if (!delta->x[0] && delta->wds <= 1) {
2180 delta = lshift(delta,Log2P);
2181 if (cmp(delta, bs) > 0)
2188 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2190 #ifdef Avoid_Underflow
2191 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2192 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2196 word0(rv) = (word0(rv) & Exp_mask)
2203 #ifdef Avoid_Underflow
2209 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2212 #ifdef Sudden_Underflow
2213 L = word0(rv) & Exp_mask;
2217 #ifdef Avoid_Underflow
2218 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2226 #ifdef Avoid_Underflow
2228 L = word0(rv) & Exp_mask;
2229 if (L <= (2*P+1)*Exp_msk1) {
2230 if (L > (P+2)*Exp_msk1)
2239 L = (word0(rv) & Exp_mask) - Exp_msk1;
2241 word0(rv) = L | Bndry_mask1;
2242 word1(rv) = 0xffffffff;
2249 #ifndef ROUND_BIASED
2250 if (!(word1(rv) & LSB))
2254 dval(rv) += ulp(dval(rv));
2255 #ifndef ROUND_BIASED
2257 dval(rv) -= ulp(dval(rv));
2258 #ifndef Sudden_Underflow
2263 #ifdef Avoid_Underflow
2269 if ((aadj = ratio(delta, bs)) <= 2.) {
2271 aadj = dval(aadj1) = 1.;
2272 else if (word1(rv) || word0(rv) & Bndry_mask) {
2273 #ifndef Sudden_Underflow
2274 if (word1(rv) == Tiny1 && !word0(rv))
2284 if (aadj < 2./FLT_RADIX)
2285 aadj = 1./FLT_RADIX;
2288 dval(aadj1) = -aadj;
2293 dval(aadj1) = dsign ? aadj : -aadj;
2294 #ifdef Check_FLT_ROUNDS
2304 if (Flt_Rounds == 0)
2308 y = word0(rv) & Exp_mask;
2312 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2313 dval(rv0) = dval(rv);
2314 word0(rv) -= P*Exp_msk1;
2315 adj = dval(aadj1) * ulp(dval(rv));
2317 if ((word0(rv) & Exp_mask) >=
2318 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2319 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2326 word0(rv) += P*Exp_msk1;
2329 #ifdef Avoid_Underflow
2330 if (scale && y <= 2*P*Exp_msk1) {
2331 if (aadj <= 0x7fffffff) {
2332 if ((z = (
int)aadj) <= 0)
2335 dval(aadj1) = dsign ? aadj : -aadj;
2337 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2339 adj = dval(aadj1) * ulp(dval(rv));
2342 #ifdef Sudden_Underflow
2343 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2344 dval(rv0) = dval(rv);
2345 word0(rv) += P*Exp_msk1;
2346 adj = dval(aadj1) * ulp(dval(rv));
2349 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2351 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2354 if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
2361 word0(rv) -= P*Exp_msk1;
2364 adj = dval(aadj1) * ulp(dval(rv));
2375 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2376 dval(aadj1) = (double)(
int)(aadj + 0.5);
2378 dval(aadj1) = -dval(aadj1);
2380 adj = dval(aadj1) * ulp(dval(rv));
2385 z = word0(rv) & Exp_mask;
2387 #ifdef Avoid_Underflow
2395 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2396 if (aadj < .4999999 || aadj > .5000001)
2399 else if (aadj < .4999999/FLT_RADIX)
2412 word0(rv0) = Exp_1 + (70 << Exp_shift);
2417 else if (!oldinexact)
2420 #ifdef Avoid_Underflow
2422 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2424 dval(rv) *= dval(rv0);
2427 if (word0(rv) == 0 && word1(rv) == 0)
2433 if (inexact && !(word0(rv) & Exp_mask)) {
2436 dval(rv0) *= dval(rv0);
2448 return sign ? -dval(rv) : dval(rv);
2451 NO_SANITIZE(
"unsigned-integer-overflow",
static int quorem(
Bigint *b,
Bigint *S));
2456 ULong *bx, *bxe, q, *sx, *sxe;
2458 ULLong borrow, carry, y, ys;
2460 ULong borrow, carry, y, ys;
2469 Bug(
"oversize b in quorem");
2477 q = *bxe / (*sxe + 1);
2480 Bug(
"oversized quotient in quorem");
2487 ys = *sx++ * (ULLong)q + carry;
2489 y = *bx - (ys & FFFFFFFF) - borrow;
2490 borrow = y >> 32 & (ULong)1;
2491 *bx++ = (ULong)(y & FFFFFFFF);
2495 ys = (si & 0xffff) * q + carry;
2496 zs = (si >> 16) * q + (ys >> 16);
2498 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2499 borrow = (y & 0x10000) >> 16;
2500 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2501 borrow = (z & 0x10000) >> 16;
2504 ys = *sx++ * q + carry;
2506 y = *bx - (ys & 0xffff) - borrow;
2507 borrow = (y & 0x10000) >> 16;
2511 }
while (sx <= sxe);
2514 while (--bxe > bx && !*bxe)
2519 if (cmp(b, S) >= 0) {
2529 y = *bx - (ys & FFFFFFFF) - borrow;
2530 borrow = y >> 32 & (ULong)1;
2531 *bx++ = (ULong)(y & FFFFFFFF);
2535 ys = (si & 0xffff) + carry;
2536 zs = (si >> 16) + (ys >> 16);
2538 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2539 borrow = (y & 0x10000) >> 16;
2540 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2541 borrow = (z & 0x10000) >> 16;
2546 y = *bx - (ys & 0xffff) - borrow;
2547 borrow = (y & 0x10000) >> 16;
2551 }
while (sx <= sxe);
2555 while (--bxe > bx && !*bxe)
2563 #ifndef MULTIPLE_THREADS
2564 static char *dtoa_result;
2567 #ifndef MULTIPLE_THREADS
2571 return dtoa_result = MALLOC(i);
2574 #define rv_alloc(i) MALLOC(i)
2578 nrv_alloc(
const char *s,
char **rve,
size_t n)
2582 t = rv = rv_alloc(n);
2583 while ((*t = *s++) != 0) t++;
2589 #define rv_strdup(s, rve) nrv_alloc((s), (rve), strlen(s)+1)
2591 #ifndef MULTIPLE_THREADS
2605 static const char INFSTR[] =
"Infinity";
2606 static const char NANSTR[] =
"NaN";
2607 static const char ZEROSTR[] =
"0";
2644 dtoa(
double d_,
int mode,
int ndigits,
int *decpt,
int *sign,
char **rve)
2680 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2681 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2682 spec_case, try_quick, half = 0;
2684 #ifndef Sudden_Underflow
2688 Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S;
2692 #ifdef Honor_FLT_ROUNDS
2696 int inexact, oldinexact;
2701 #ifndef MULTIPLE_THREADS
2703 freedtoa(dtoa_result);
2708 if (word0(d) & Sign_bit) {
2711 word0(d) &= ~Sign_bit;
2716 #if defined(IEEE_Arith) + defined(VAX)
2718 if ((word0(d) & Exp_mask) == Exp_mask)
2720 if (word0(d) == 0x8000)
2726 if (!word1(d) && !(word0(d) & 0xfffff))
2727 return rv_strdup(INFSTR, rve);
2729 return rv_strdup(NANSTR, rve);
2737 return rv_strdup(ZEROSTR, rve);
2741 try_quick = oldinexact = get_inexact();
2744 #ifdef Honor_FLT_ROUNDS
2745 if ((rounding = Flt_Rounds) >= 2) {
2747 rounding = rounding == 2 ? 0 : 2;
2754 b = d2b(dval(d), &be, &bbits);
2755 #ifdef Sudden_Underflow
2756 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2758 if ((i = (
int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
2761 word0(d2) &= Frac_mask1;
2762 word0(d2) |= Exp_11;
2764 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2795 #ifndef Sudden_Underflow
2801 i = bbits + be + (Bias + (P-1) - 1);
2802 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2803 : word1(d) << (32 - i);
2805 word0(d2) -= 31*Exp_msk1;
2806 i -= (Bias + (P-1) - 1) + 1;
2810 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2812 if (ds < 0. && ds != k)
2815 if (k >= 0 && k <= Ten_pmax) {
2816 if (dval(d) < tens[k])
2839 if (mode < 0 || mode > 9)
2843 #ifdef Check_FLT_ROUNDS
2844 try_quick = Rounding == 1;
2868 ilim = ilim1 = i = ndigits;
2874 i = ndigits + k + 1;
2880 s = s0 = rv_alloc(i+1);
2882 #ifdef Honor_FLT_ROUNDS
2883 if (mode > 1 && rounding != 1)
2887 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2902 dval(d) /= bigtens[n_bigtens-1];
2905 for (; j; j >>= 1, i++)
2912 else if ((j1 = -k) != 0) {
2913 dval(d) *= tens[j1 & 0xf];
2914 for (j = j1 >> 4; j; j >>= 1, i++)
2917 dval(d) *= bigtens[i];
2920 if (k_check && dval(d) < 1. && ilim > 0) {
2928 dval(eps) = ieps*dval(d) + 7.;
2929 word0(eps) -= (P-1)*Exp_msk1;
2933 if (dval(d) > dval(eps))
2935 if (dval(d) < -dval(eps))
2939 #ifndef No_leftright
2944 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2948 *s++ =
'0' + (int)L;
2949 if (dval(d) < dval(eps))
2951 if (1. - dval(d) < dval(eps))
2962 dval(eps) *= tens[ilim-1];
2963 for (i = 1;; i++, dval(d) *= 10.) {
2964 L = (Long)(dval(d));
2965 if (!(dval(d) -= L))
2967 *s++ =
'0' + (int)L;
2969 if (dval(d) > 0.5 + dval(eps))
2971 else if (dval(d) < 0.5 - dval(eps)) {
2972 while (*--s ==
'0') ;
2977 if ((*(s-1) -
'0') & 1) {
2983 #ifndef No_leftright
2995 if (be >= 0 && k <= Int_max) {
2998 if (ndigits < 0 && ilim <= 0) {
3000 if (ilim < 0 || dval(d) <= 5*ds)
3004 for (i = 1;; i++, dval(d) *= 10.) {
3005 L = (Long)(dval(d) / ds);
3007 #ifdef Check_FLT_ROUNDS
3014 *s++ =
'0' + (int)L;
3022 #ifdef Honor_FLT_ROUNDS
3026 case 2:
goto bump_up;
3030 if (dval(d) > ds || (dval(d) == ds && (L & 1))) {
3050 #ifndef Sudden_Underflow
3051 denorm ? be + (Bias + (P-1) - 1 + 1) :
3054 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3062 if (m2 > 0 && s2 > 0) {
3063 i = m2 < s2 ? m2 : s2;
3071 mhi = pow5mult(mhi, m5);
3076 if ((j = b5 - m5) != 0)
3080 b = pow5mult(b, b5);
3084 S = pow5mult(S, s5);
3089 if ((mode < 2 || leftright)
3090 #ifdef Honor_FLT_ROUNDS
3094 if (!word1(d) && !(word0(d) & Bndry_mask)
3095 #ifndef Sudden_Underflow
3096 && word0(d) & (Exp_mask & ~Exp_msk1)
3114 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
3117 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0)
3139 b = multadd(b, 10, 0);
3141 mhi = multadd(mhi, 10, 0);
3145 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3146 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3159 mhi = lshift(mhi, m2);
3167 mhi = Balloc(mhi->k);
3169 mhi = lshift(mhi, Log2P);
3173 dig = quorem(b,S) +
'0';
3178 delta = diff(S, mhi);
3179 j1 = delta->sign ? 1 : cmp(b, delta);
3181 #ifndef ROUND_BIASED
3182 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3183 #ifdef Honor_FLT_ROUNDS
3192 else if (!b->x[0] && b->wds <= 1)
3199 if (j < 0 || (j == 0 && mode != 1
3200 #ifndef ROUND_BIASED
3204 if (!b->x[0] && b->wds <= 1) {
3210 #ifdef Honor_FLT_ROUNDS
3213 case 0:
goto accept_dig;
3214 case 2:
goto keep_dig;
3220 if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ ==
'9')
3228 #ifdef Honor_FLT_ROUNDS
3240 #ifdef Honor_FLT_ROUNDS
3246 b = multadd(b, 10, 0);
3248 mlo = mhi = multadd(mhi, 10, 0);
3250 mlo = multadd(mlo, 10, 0);
3251 mhi = multadd(mhi, 10, 0);
3257 *s++ = dig = quorem(b,S) +
'0';
3258 if (!b->x[0] && b->wds <= 1) {
3266 b = multadd(b, 10, 0);
3271 #ifdef Honor_FLT_ROUNDS
3273 case 0:
goto trimzeros;
3274 case 2:
goto roundoff;
3279 if (j > 0 || (j == 0 && (dig & 1))) {
3287 if (!half || (*s -
'0') & 1)
3291 while (*--s ==
'0') ;
3297 if (mlo && mlo != mhi)
3305 word0(d) = Exp_1 + (70 << Exp_shift);
3310 else if (!oldinexact)
3347 #define DBL_MANH_SIZE 20
3348 #define DBL_MANL_SIZE 32
3349 #define DBL_ADJ (DBL_MAX_EXP - 2)
3350 #define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
3351 #define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1)
3352 #define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | ((v) << Exp_shift)))
3353 #define dmanh_get(u) ((uint32_t)(word0(u) & Frac_mask))
3354 #define dmanl_get(u) ((uint32_t)word1(u))
3382 hdtoa(
double d,
const char *xdigs,
int ndigits,
int *decpt,
int *sign,
char **rve)
3387 uint32_t manh, manl;
3390 if (word0(u) & Sign_bit) {
3393 word0(u) &= ~Sign_bit;
3400 return rv_strdup(INFSTR, rve);
3402 else if (isnan(d)) {
3404 return rv_strdup(NANSTR, rve);
3406 else if (d == 0.0) {
3408 return rv_strdup(ZEROSTR, rve);
3410 else if (dexp_get(u)) {
3411 *decpt = dexp_get(u) - DBL_ADJ;
3414 u.d *= 5.363123171977039e+154 ;
3415 *decpt = dexp_get(u) - (514 + DBL_ADJ);
3425 bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
3426 s0 = rv_alloc(bufsize+1);
3429 if (SIGFIGS > ndigits && ndigits > 0) {
3431 int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
3432 dexp_set(u, offset);
3435 *decpt += dexp_get(u) - offset;
3438 manh = dmanh_get(u);
3439 manl = dmanl_get(u);
3441 for (s = s0 + 1; s < s0 + bufsize; s++) {
3442 *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
3443 manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
3449 for (ndigits = SIGFIGS; s0[ndigits - 1] ==
'0'; ndigits--)
#define ISDIGIT
Old name of rb_isdigit.
#define ASSUME
Old name of RBIMPL_ASSUME.
int len
Length of the buffer.
#define strtod(s, e)
Just another name of ruby_strtod.
#define errno
Ractor-aware version of errno.