OLD | NEW |
1 /**************************************************************** | 1 /**************************************************************** |
2 * | 2 * |
3 * The author of this software is David M. Gay. | 3 * The author of this software is David M. Gay. |
4 * | 4 * |
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. | 5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. |
6 * | 6 * |
7 * Permission to use, copy, modify, and distribute this software for any | 7 * Permission to use, copy, modify, and distribute this software for any |
8 * purpose without fee is hereby granted, provided that this entire notice | 8 * purpose without fee is hereby granted, provided that this entire notice |
9 * is included in all copies of any software which is or includes a copy | 9 * is included in all copies of any software which is or includes a copy |
10 * or modification of this software and in all copies of the supporting | 10 * or modification of this software and in all copies of the supporting |
(...skipping 14 matching lines...) Expand all Loading... |
25 * before invoking strtod or dtoa. If the machine uses (the equivalent | 25 * before invoking strtod or dtoa. If the machine uses (the equivalent |
26 * of) Intel 80x87 arithmetic, the call | 26 * of) Intel 80x87 arithmetic, the call |
27 * _control87(PC_53, MCW_PC); | 27 * _control87(PC_53, MCW_PC); |
28 * does this with many compilers. Whether this or another call is | 28 * does this with many compilers. Whether this or another call is |
29 * appropriate depends on the compiler; for this to work, it may be | 29 * appropriate depends on the compiler; for this to work, it may be |
30 * necessary to #include "float.h" or another system-dependent header | 30 * necessary to #include "float.h" or another system-dependent header |
31 * file. | 31 * file. |
32 */ | 32 */ |
33 | 33 |
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. | 34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. |
| 35 * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.) |
35 * | 36 * |
36 * This strtod returns a nearest machine number to the input decimal | 37 * This strtod returns a nearest machine number to the input decimal |
37 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are | 38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are |
38 * broken by the IEEE round-even rule. Otherwise ties are broken by | 39 * broken by the IEEE round-even rule. Otherwise ties are broken by |
39 * biased rounding (add half and chop). | 40 * biased rounding (add half and chop). |
40 * | 41 * |
41 * Inspired loosely by William D. Clinger's paper "How to Read Floating | 42 * Inspired loosely by William D. Clinger's paper "How to Read Floating |
42 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. | 43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. |
43 * | 44 * |
44 * Modifications: | 45 * Modifications: |
(...skipping 18 matching lines...) Expand all Loading... |
63 | 64 |
64 /* | 65 /* |
65 * #define IEEE_8087 for IEEE-arithmetic machines where the least | 66 * #define IEEE_8087 for IEEE-arithmetic machines where the least |
66 * significant byte has the lowest address. | 67 * significant byte has the lowest address. |
67 * #define IEEE_MC68k for IEEE-arithmetic machines where the most | 68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most |
68 * significant byte has the lowest address. | 69 * significant byte has the lowest address. |
69 * #define Long int on machines with 32-bit ints and 64-bit longs. | 70 * #define Long int on machines with 32-bit ints and 64-bit longs. |
70 * #define IBM for IBM mainframe-style floating-point arithmetic. | 71 * #define IBM for IBM mainframe-style floating-point arithmetic. |
71 * #define VAX for VAX-style floating-point arithmetic (D_floating). | 72 * #define VAX for VAX-style floating-point arithmetic (D_floating). |
72 * #define No_leftright to omit left-right logic in fast floating-point | 73 * #define No_leftright to omit left-right logic in fast floating-point |
73 *» computation of dtoa. | 74 *» computation of dtoa. This will cause dtoa modes 4 and 5 to be |
| 75 *» treated the same as modes 2 and 3 for some inputs. |
74 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 | 76 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 |
75 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS | 77 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS |
76 * is also #defined, fegetround() will be queried for the rounding mode. | 78 * is also #defined, fegetround() will be queried for the rounding mode. |
77 * Note that both FLT_ROUNDS and fegetround() are specified by the C99 | 79 * Note that both FLT_ROUNDS and fegetround() are specified by the C99 |
78 * standard (and are specified to be consistent, with fesetround() | 80 * standard (and are specified to be consistent, with fesetround() |
79 * affecting the value of FLT_ROUNDS), but that some (Linux) systems | 81 * affecting the value of FLT_ROUNDS), but that some (Linux) systems |
80 * do not work correctly in this regard, so using fegetround() is more | 82 * do not work correctly in this regard, so using fegetround() is more |
81 *» portable than using FLT_FOUNDS directly. | 83 *» portable than using FLT_ROUNDS directly. |
82 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 | 84 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 |
83 * and Honor_FLT_ROUNDS is not #defined. | 85 * and Honor_FLT_ROUNDS is not #defined. |
84 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines | 86 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines |
85 * that use extended-precision instructions to compute rounded | 87 * that use extended-precision instructions to compute rounded |
86 * products and quotients) with IBM. | 88 * products and quotients) with IBM. |
87 * #define ROUND_BIASED for IEEE-format with biased rounding. | 89 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic |
| 90 *» that rounds toward +Infinity. |
| 91 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased |
| 92 *» rounding when the underlying floating-point arithmetic uses |
| 93 *» unbiased rounding. This prevent using ordinary floating-point |
| 94 *» arithmetic when the result could be computed with one rounding error. |
88 * #define Inaccurate_Divide for IEEE-format with correctly rounded | 95 * #define Inaccurate_Divide for IEEE-format with correctly rounded |
89 * products but inaccurate quotients, e.g., for Intel i860. | 96 * products but inaccurate quotients, e.g., for Intel i860. |
90 * #define NO_LONG_LONG on machines that do not have a "long long" | 97 * #define NO_LONG_LONG on machines that do not have a "long long" |
91 * integer type (of >= 64 bits). On such machines, you can | 98 * integer type (of >= 64 bits). On such machines, you can |
92 * #define Just_16 to store 16 bits per 32-bit Long when doing | 99 * #define Just_16 to store 16 bits per 32-bit Long when doing |
93 * high-precision integer arithmetic. Whether this speeds things | 100 * high-precision integer arithmetic. Whether this speeds things |
94 * up or slows things down depends on the machine and the number | 101 * up or slows things down depends on the machine and the number |
95 * being converted. If long long is available and the name is | 102 * being converted. If long long is available and the name is |
96 * something other than "long long", #define Llong to be the name, | 103 * something other than "long long", #define Llong to be the name, |
97 * and if "unsigned Llong" does not work as an unsigned version of | 104 * and if "unsigned Llong" does not work as an unsigned version of |
(...skipping 353 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
451 #define Log2P 1 | 458 #define Log2P 1 |
452 #define Tiny0 0x80 | 459 #define Tiny0 0x80 |
453 #define Tiny1 0 | 460 #define Tiny1 0 |
454 #define Quick_max 15 | 461 #define Quick_max 15 |
455 #define Int_max 15 | 462 #define Int_max 15 |
456 #endif /* IBM, VAX */ | 463 #endif /* IBM, VAX */ |
457 #endif /* IEEE_Arith */ | 464 #endif /* IEEE_Arith */ |
458 | 465 |
459 #ifndef IEEE_Arith | 466 #ifndef IEEE_Arith |
460 #define ROUND_BIASED | 467 #define ROUND_BIASED |
| 468 #else |
| 469 #ifdef ROUND_BIASED_without_Round_Up |
| 470 #undef ROUND_BIASED |
| 471 #define ROUND_BIASED |
| 472 #endif |
461 #endif | 473 #endif |
462 | 474 |
463 #ifdef RND_PRODQUOT | 475 #ifdef RND_PRODQUOT |
464 #define rounded_product(a,b) a = rnd_prod(a, b) | 476 #define rounded_product(a,b) a = rnd_prod(a, b) |
465 #define rounded_quotient(a,b) a = rnd_quot(a, b) | 477 #define rounded_quotient(a,b) a = rnd_quot(a, b) |
466 #ifdef KR_headers | 478 #ifdef KR_headers |
467 extern double rnd_prod(), rnd_quot(); | 479 extern double rnd_prod(), rnd_quot(); |
468 #else | 480 #else |
469 extern double rnd_prod(double, double), rnd_quot(double, double); | 481 extern double rnd_prod(double, double), rnd_quot(double, double); |
470 #endif | 482 #endif |
(...skipping 185 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
656 b->wds = wds; | 668 b->wds = wds; |
657 } | 669 } |
658 return b; | 670 return b; |
659 } | 671 } |
660 | 672 |
661 static Bigint * | 673 static Bigint * |
662 s2b | 674 s2b |
663 #ifdef KR_headers | 675 #ifdef KR_headers |
664 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9; | 676 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9; |
665 #else | 677 #else |
666 » (CONST char *s, int nd0, int nd, ULong y9, int dplen) | 678 » (const char *s, int nd0, int nd, ULong y9, int dplen) |
667 #endif | 679 #endif |
668 { | 680 { |
669 Bigint *b; | 681 Bigint *b; |
670 int i, k; | 682 int i, k; |
671 Long x, y; | 683 Long x, y; |
672 | 684 |
673 x = (nd + 8) / 9; | 685 x = (nd + 8) / 9; |
674 for(k = 0, y = 1; x > y; y <<= 1, k++) ; | 686 for(k = 0, y = 1; x > y; y <<= 1, k++) ; |
675 #ifdef Pack_32 | 687 #ifdef Pack_32 |
676 b = Balloc(k); | 688 b = Balloc(k); |
(...skipping 816 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1493 #endif | 1505 #endif |
1494 #endif | 1506 #endif |
1495 | 1507 |
1496 #ifndef Need_Hexdig | 1508 #ifndef Need_Hexdig |
1497 #ifndef NO_HEX_FP | 1509 #ifndef NO_HEX_FP |
1498 #define Need_Hexdig | 1510 #define Need_Hexdig |
1499 #endif | 1511 #endif |
1500 #endif | 1512 #endif |
1501 | 1513 |
1502 #ifdef Need_Hexdig /*{*/ | 1514 #ifdef Need_Hexdig /*{*/ |
| 1515 #if 0 |
1503 static unsigned char hexdig[256]; | 1516 static unsigned char hexdig[256]; |
1504 | 1517 |
1505 static void | 1518 static void |
1506 #ifdef KR_headers | |
1507 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc; | |
1508 #else | |
1509 htinit(unsigned char *h, unsigned char *s, int inc) | 1519 htinit(unsigned char *h, unsigned char *s, int inc) |
1510 #endif | |
1511 { | 1520 { |
1512 int i, j; | 1521 int i, j; |
1513 for(i = 0; (j = s[i]) !=0; i++) | 1522 for(i = 0; (j = s[i]) !=0; i++) |
1514 h[j] = (unsigned char)(i + inc); | 1523 h[j] = (unsigned char)(i + inc); |
1515 } | 1524 } |
1516 | 1525 |
1517 static void | 1526 static void |
1518 #ifdef KR_headers | 1527 hexdig_init(void)» /* Use of hexdig_init omitted 20121220 to avoid a */ |
1519 hexdig_init() | 1528 » » » /* race condition when multiple threads are used. */ |
1520 #else | |
1521 hexdig_init(void) | |
1522 #endif | |
1523 { | 1529 { |
1524 #define USC (unsigned char *) | 1530 #define USC (unsigned char *) |
1525 htinit(hexdig, USC "0123456789", 0x10); | 1531 htinit(hexdig, USC "0123456789", 0x10); |
1526 htinit(hexdig, USC "abcdef", 0x10 + 10); | 1532 htinit(hexdig, USC "abcdef", 0x10 + 10); |
1527 htinit(hexdig, USC "ABCDEF", 0x10 + 10); | 1533 htinit(hexdig, USC "ABCDEF", 0x10 + 10); |
1528 } | 1534 } |
| 1535 #else |
| 1536 static unsigned char hexdig[256] = { |
| 1537 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1538 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1539 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1540 16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0, |
| 1541 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, |
| 1542 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1543 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, |
| 1544 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1545 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1546 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1547 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1548 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1549 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1550 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1551 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, |
| 1552 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 |
| 1553 }; |
| 1554 #endif |
1529 #endif /* } Need_Hexdig */ | 1555 #endif /* } Need_Hexdig */ |
1530 | 1556 |
1531 #ifdef INFNAN_CHECK | 1557 #ifdef INFNAN_CHECK |
1532 | 1558 |
1533 #ifndef NAN_WORD0 | 1559 #ifndef NAN_WORD0 |
1534 #define NAN_WORD0 0x7ff80000 | 1560 #define NAN_WORD0 0x7ff80000 |
1535 #endif | 1561 #endif |
1536 | 1562 |
1537 #ifndef NAN_WORD1 | 1563 #ifndef NAN_WORD1 |
1538 #define NAN_WORD1 0 | 1564 #define NAN_WORD1 0 |
1539 #endif | 1565 #endif |
1540 | 1566 |
1541 static int | 1567 static int |
1542 match | 1568 match |
1543 #ifdef KR_headers | 1569 #ifdef KR_headers |
1544 (sp, t) char **sp, *t; | 1570 (sp, t) char **sp, *t; |
1545 #else | 1571 #else |
1546 » (CONST char **sp, CONST char *t) | 1572 » (const char **sp, const char *t) |
1547 #endif | 1573 #endif |
1548 { | 1574 { |
1549 int c, d; | 1575 int c, d; |
1550 CONST char *s = *sp; | 1576 CONST char *s = *sp; |
1551 | 1577 |
1552 for(d = *t++; d; d = *t++) { | 1578 for(d = *t++; d; d = *t++) { |
1553 if ((c = *++s) >= 'A' && c <= 'Z') | 1579 if ((c = *++s) >= 'A' && c <= 'Z') |
1554 c += 'a' - 'A'; | 1580 c += 'a' - 'A'; |
1555 if (c != d) | 1581 if (c != d) |
1556 return 0; | 1582 return 0; |
1557 } | 1583 } |
1558 *sp = s + 1; | 1584 *sp = s + 1; |
1559 return 1; | 1585 return 1; |
1560 } | 1586 } |
1561 | 1587 |
1562 #ifndef No_Hex_NaN | 1588 #ifndef No_Hex_NaN |
1563 static void | 1589 static void |
1564 hexnan | 1590 hexnan |
1565 #ifdef KR_headers | 1591 #ifdef KR_headers |
1566 (rvp, sp) U *rvp; CONST char **sp; | 1592 (rvp, sp) U *rvp; CONST char **sp; |
1567 #else | 1593 #else |
1568 » (U *rvp, CONST char **sp) | 1594 » (U *rvp, const char **sp) |
1569 #endif | 1595 #endif |
1570 { | 1596 { |
1571 ULong c, x[2]; | 1597 ULong c, x[2]; |
1572 CONST char *s; | 1598 CONST char *s; |
1573 int c1, havedig, udx0, xshift; | 1599 int c1, havedig, udx0, xshift; |
1574 | 1600 |
1575 » if (!hexdig['0']) | 1601 » /**** if (!hexdig['0']) hexdig_init(); ****/ |
1576 » » hexdig_init(); | |
1577 x[0] = x[1] = 0; | 1602 x[0] = x[1] = 0; |
1578 havedig = xshift = 0; | 1603 havedig = xshift = 0; |
1579 udx0 = 1; | 1604 udx0 = 1; |
1580 s = *sp; | 1605 s = *sp; |
1581 /* allow optional initial 0x or 0X */ | 1606 /* allow optional initial 0x or 0X */ |
1582 for(c = *(CONST unsigned char*)(s+1); c && c <= ' '; c = *(CONST unsigne
d char*)(s+1)) | 1607 for(c = *(CONST unsigned char*)(s+1); c && c <= ' '; c = *(CONST unsigne
d char*)(s+1)) |
1583 ++s; | 1608 ++s; |
1584 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) | 1609 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')) |
1585 s += 2; | 1610 s += 2; |
1586 for(c = *(CONST unsigned char*)++s; c; c = *(CONST unsigned char*)++s) { | 1611 for(c = *(CONST unsigned char*)++s; c; c = *(CONST unsigned char*)++s) { |
(...skipping 46 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1633 | 1658 |
1634 #ifdef Pack_32 | 1659 #ifdef Pack_32 |
1635 #define ULbits 32 | 1660 #define ULbits 32 |
1636 #define kshift 5 | 1661 #define kshift 5 |
1637 #define kmask 31 | 1662 #define kmask 31 |
1638 #else | 1663 #else |
1639 #define ULbits 16 | 1664 #define ULbits 16 |
1640 #define kshift 4 | 1665 #define kshift 4 |
1641 #define kmask 15 | 1666 #define kmask 15 |
1642 #endif | 1667 #endif |
| 1668 |
| 1669 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/ |
| 1670 static Bigint * |
| 1671 #ifdef KR_headers |
| 1672 increment(b) Bigint *b; |
| 1673 #else |
| 1674 increment(Bigint *b) |
| 1675 #endif |
| 1676 { |
| 1677 ULong *x, *xe; |
| 1678 Bigint *b1; |
| 1679 |
| 1680 x = b->x; |
| 1681 xe = x + b->wds; |
| 1682 do { |
| 1683 if (*x < (ULong)0xffffffffL) { |
| 1684 ++*x; |
| 1685 return b; |
| 1686 } |
| 1687 *x++ = 0; |
| 1688 } while(x < xe); |
| 1689 { |
| 1690 if (b->wds >= b->maxwds) { |
| 1691 b1 = Balloc(b->k+1); |
| 1692 Bcopy(b1,b); |
| 1693 Bfree(b); |
| 1694 b = b1; |
| 1695 } |
| 1696 b->x[b->wds++] = 1; |
| 1697 } |
| 1698 return b; |
| 1699 } |
| 1700 |
| 1701 #endif /*}*/ |
| 1702 |
1643 #ifndef NO_HEX_FP /*{*/ | 1703 #ifndef NO_HEX_FP /*{*/ |
1644 | 1704 |
1645 static void | 1705 static void |
1646 #ifdef KR_headers | 1706 #ifdef KR_headers |
1647 rshift(b, k) Bigint *b; int k; | 1707 rshift(b, k) Bigint *b; int k; |
1648 #else | 1708 #else |
1649 rshift(Bigint *b, int k) | 1709 rshift(Bigint *b, int k) |
1650 #endif | 1710 #endif |
1651 { | 1711 { |
1652 ULong *x, *x1, *xe, y; | 1712 ULong *x, *x1, *xe, y; |
(...skipping 52 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1705 return 0; | 1765 return 0; |
1706 } | 1766 } |
1707 | 1767 |
1708 enum { /* rounding values: same as FLT_ROUNDS */ | 1768 enum { /* rounding values: same as FLT_ROUNDS */ |
1709 Round_zero = 0, | 1769 Round_zero = 0, |
1710 Round_near = 1, | 1770 Round_near = 1, |
1711 Round_up = 2, | 1771 Round_up = 2, |
1712 Round_down = 3 | 1772 Round_down = 3 |
1713 }; | 1773 }; |
1714 | 1774 |
1715 static Bigint * | |
1716 #ifdef KR_headers | |
1717 increment(b) Bigint *b; | |
1718 #else | |
1719 increment(Bigint *b) | |
1720 #endif | |
1721 { | |
1722 ULong *x, *xe; | |
1723 Bigint *b1; | |
1724 | |
1725 x = b->x; | |
1726 xe = x + b->wds; | |
1727 do { | |
1728 if (*x < (ULong)0xffffffffL) { | |
1729 ++*x; | |
1730 return b; | |
1731 } | |
1732 *x++ = 0; | |
1733 } while(x < xe); | |
1734 { | |
1735 if (b->wds >= b->maxwds) { | |
1736 b1 = Balloc(b->k+1); | |
1737 Bcopy(b1,b); | |
1738 Bfree(b); | |
1739 b = b1; | |
1740 } | |
1741 b->x[b->wds++] = 1; | |
1742 } | |
1743 return b; | |
1744 } | |
1745 | |
1746 void | 1775 void |
1747 #ifdef KR_headers | 1776 #ifdef KR_headers |
1748 gethex(sp, rvp, rounding, sign) | 1777 gethex(sp, rvp, rounding, sign) |
1749 CONST char **sp; U *rvp; int rounding, sign; | 1778 CONST char **sp; U *rvp; int rounding, sign; |
1750 #else | 1779 #else |
1751 gethex( CONST char **sp, U *rvp, int rounding, int sign) | 1780 gethex( CONST char **sp, U *rvp, int rounding, int sign) |
1752 #endif | 1781 #endif |
1753 { | 1782 { |
1754 Bigint *b; | 1783 Bigint *b; |
1755 CONST unsigned char *decpt, *s0, *s, *s1; | 1784 CONST unsigned char *decpt, *s0, *s, *s1; |
(...skipping 30 matching lines...) Expand all Loading... |
1786 if ((decimalpoint_cache = (unsigned char*) | 1815 if ((decimalpoint_cache = (unsigned char*) |
1787 MALLOC(strlen((CONST char*)s0) + 1))) { | 1816 MALLOC(strlen((CONST char*)s0) + 1))) { |
1788 strcpy((char*)decimalpoint_cache, (CONST char*)s0); | 1817 strcpy((char*)decimalpoint_cache, (CONST char*)s0); |
1789 s0 = decimalpoint_cache; | 1818 s0 = decimalpoint_cache; |
1790 } | 1819 } |
1791 } | 1820 } |
1792 decimalpoint = s0; | 1821 decimalpoint = s0; |
1793 #endif | 1822 #endif |
1794 #endif | 1823 #endif |
1795 | 1824 |
1796 » if (!hexdig['0']) | 1825 » /**** if (!hexdig['0']) hexdig_init(); ****/ |
1797 » » hexdig_init(); | |
1798 havedig = 0; | 1826 havedig = 0; |
1799 s0 = *(CONST unsigned char **)sp + 2; | 1827 s0 = *(CONST unsigned char **)sp + 2; |
1800 while(s0[havedig] == '0') | 1828 while(s0[havedig] == '0') |
1801 havedig++; | 1829 havedig++; |
1802 s0 += havedig; | 1830 s0 += havedig; |
1803 s = s0; | 1831 s = s0; |
1804 decpt = 0; | 1832 decpt = 0; |
1805 zret = 0; | 1833 zret = 0; |
1806 e = 0; | 1834 e = 0; |
1807 if (hexdig[*s]) | 1835 if (hexdig[*s]) |
(...skipping 79 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1887 break; | 1915 break; |
1888 goto ret_tiny; | 1916 goto ret_tiny; |
1889 case Round_down: | 1917 case Round_down: |
1890 if (!sign) | 1918 if (!sign) |
1891 break; | 1919 break; |
1892 goto ret_tiny; | 1920 goto ret_tiny; |
1893 } | 1921 } |
1894 #endif | 1922 #endif |
1895 goto retz; | 1923 goto retz; |
1896 #ifdef IEEE_Arith | 1924 #ifdef IEEE_Arith |
| 1925 ret_tinyf: |
| 1926 Bfree(b); |
1897 ret_tiny: | 1927 ret_tiny: |
1898 #ifndef NO_ERRNO | 1928 #ifndef NO_ERRNO |
1899 errno = ERANGE; | 1929 errno = ERANGE; |
1900 #endif | 1930 #endif |
1901 word0(rvp) = 0; | 1931 word0(rvp) = 0; |
1902 word1(rvp) = 1; | 1932 word1(rvp) = 1; |
1903 return; | 1933 return; |
1904 #endif /* IEEE_Arith */ | 1934 #endif /* IEEE_Arith */ |
1905 } | 1935 } |
1906 switch(rounding) { | 1936 switch(rounding) { |
(...skipping 80 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1987 } | 2017 } |
1988 denorm = 0; | 2018 denorm = 0; |
1989 if (e < emin) { | 2019 if (e < emin) { |
1990 denorm = 1; | 2020 denorm = 1; |
1991 n = emin - e; | 2021 n = emin - e; |
1992 if (n >= nbits) { | 2022 if (n >= nbits) { |
1993 #ifdef IEEE_Arith /*{*/ | 2023 #ifdef IEEE_Arith /*{*/ |
1994 switch (rounding) { | 2024 switch (rounding) { |
1995 case Round_near: | 2025 case Round_near: |
1996 if (n == nbits && (n < 2 || any_on(b,n-1))) | 2026 if (n == nbits && (n < 2 || any_on(b,n-1))) |
1997 » » » » » goto ret_tiny; | 2027 » » » » » goto ret_tinyf; |
1998 break; | 2028 break; |
1999 case Round_up: | 2029 case Round_up: |
2000 if (!sign) | 2030 if (!sign) |
2001 » » » » » goto ret_tiny; | 2031 » » » » » goto ret_tinyf; |
2002 break; | 2032 break; |
2003 case Round_down: | 2033 case Round_down: |
2004 if (sign) | 2034 if (sign) |
2005 » » » » » goto ret_tiny; | 2035 » » » » » goto ret_tinyf; |
2006 } | 2036 } |
2007 #endif /* } IEEE_Arith */ | 2037 #endif /* } IEEE_Arith */ |
2008 Bfree(b); | 2038 Bfree(b); |
2009 retz: | 2039 retz: |
2010 #ifndef NO_ERRNO | 2040 #ifndef NO_ERRNO |
2011 errno = ERANGE; | 2041 errno = ERANGE; |
2012 #endif | 2042 #endif |
2013 retz1: | 2043 retz1: |
2014 rvp->d = 0.; | 2044 rvp->d = 0.; |
2015 return; | 2045 return; |
(...skipping 79 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2095 #endif | 2125 #endif |
2096 #ifdef VAX | 2126 #ifdef VAX |
2097 /* The next two lines ignore swap of low- and high-order 2 bytes. */ | 2127 /* The next two lines ignore swap of low- and high-order 2 bytes. */ |
2098 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */ | 2128 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */ |
2099 /* word1(rvp) = b->x[0]; */ | 2129 /* word1(rvp) = b->x[0]; */ |
2100 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b-
>x[1] << 16); | 2130 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b-
>x[1] << 16); |
2101 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16); | 2131 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16); |
2102 #endif | 2132 #endif |
2103 Bfree(b); | 2133 Bfree(b); |
2104 } | 2134 } |
2105 #endif /*}!NO_HEX_FP*/ | 2135 #endif /*!NO_HEX_FP}*/ |
2106 | 2136 |
2107 static int | 2137 static int |
2108 #ifdef KR_headers | 2138 #ifdef KR_headers |
2109 dshift(b, p2) Bigint *b; int p2; | 2139 dshift(b, p2) Bigint *b; int p2; |
2110 #else | 2140 #else |
2111 dshift(Bigint *b, int p2) | 2141 dshift(Bigint *b, int p2) |
2112 #endif | 2142 #endif |
2113 { | 2143 { |
2114 int rv = hi0bits(b->x[b->wds-1]) - 4; | 2144 int rv = hi0bits(b->x[b->wds-1]) - 4; |
2115 if (p2 > 0) | 2145 if (p2 > 0) |
(...skipping 26 matching lines...) Expand all Loading... |
2142 /*debug*/ Bug("oversize b in quorem"); | 2172 /*debug*/ Bug("oversize b in quorem"); |
2143 #endif | 2173 #endif |
2144 if (b->wds < n) | 2174 if (b->wds < n) |
2145 return 0; | 2175 return 0; |
2146 sx = S->x; | 2176 sx = S->x; |
2147 sxe = sx + --n; | 2177 sxe = sx + --n; |
2148 bx = b->x; | 2178 bx = b->x; |
2149 bxe = bx + n; | 2179 bxe = bx + n; |
2150 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ | 2180 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ |
2151 #ifdef DEBUG | 2181 #ifdef DEBUG |
| 2182 #ifdef NO_STRTOD_BIGCOMP |
2152 /*debug*/ if (q > 9) | 2183 /*debug*/ if (q > 9) |
| 2184 #else |
| 2185 /* An oversized q is possible when quorem is called from bigcomp and */ |
| 2186 /* the input is near, e.g., twice the smallest denormalized number. */ |
| 2187 /*debug*/ if (q > 15) |
| 2188 #endif |
2153 /*debug*/ Bug("oversized quotient in quorem"); | 2189 /*debug*/ Bug("oversized quotient in quorem"); |
2154 #endif | 2190 #endif |
2155 if (q) { | 2191 if (q) { |
2156 borrow = 0; | 2192 borrow = 0; |
2157 carry = 0; | 2193 carry = 0; |
2158 do { | 2194 do { |
2159 #ifdef ULLong | 2195 #ifdef ULLong |
2160 ys = *sx++ * (ULLong)q + carry; | 2196 ys = *sx++ * (ULLong)q + carry; |
2161 carry = ys >> 32; | 2197 carry = ys >> 32; |
2162 y = *bx - (ys & FFFFFFFF) - borrow; | 2198 y = *bx - (ys & FFFFFFFF) - borrow; |
(...skipping 65 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2228 bxe = bx + n; | 2264 bxe = bx + n; |
2229 if (!*bxe) { | 2265 if (!*bxe) { |
2230 while(--bxe > bx && !*bxe) | 2266 while(--bxe > bx && !*bxe) |
2231 --n; | 2267 --n; |
2232 b->wds = n; | 2268 b->wds = n; |
2233 } | 2269 } |
2234 } | 2270 } |
2235 return q; | 2271 return q; |
2236 } | 2272 } |
2237 | 2273 |
| 2274 #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/ |
| 2275 static double |
| 2276 sulp |
| 2277 #ifdef KR_headers |
| 2278 (x, bc) U *x; BCinfo *bc; |
| 2279 #else |
| 2280 (U *x, BCinfo *bc) |
| 2281 #endif |
| 2282 { |
| 2283 U u; |
| 2284 double rv; |
| 2285 int i; |
| 2286 |
| 2287 rv = ulp(x); |
| 2288 if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <
= 0) |
| 2289 return rv; /* Is there an example where i <= 0 ? */ |
| 2290 word0(&u) = Exp_1 + (i << Exp_shift); |
| 2291 word1(&u) = 0; |
| 2292 return rv * u.d; |
| 2293 } |
| 2294 #endif /*}*/ |
| 2295 |
2238 #ifndef NO_STRTOD_BIGCOMP | 2296 #ifndef NO_STRTOD_BIGCOMP |
2239 | |
2240 static void | 2297 static void |
2241 bigcomp | 2298 bigcomp |
2242 #ifdef KR_headers | 2299 #ifdef KR_headers |
2243 (rv, s0, bc) | 2300 (rv, s0, bc) |
2244 U *rv; CONST char *s0; BCinfo *bc; | 2301 U *rv; CONST char *s0; BCinfo *bc; |
2245 #else | 2302 #else |
2246 » (U *rv, CONST char *s0, BCinfo *bc) | 2303 » (U *rv, const char *s0, BCinfo *bc) |
2247 #endif | 2304 #endif |
2248 { | 2305 { |
2249 Bigint *b, *d; | 2306 Bigint *b, *d; |
2250 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; | 2307 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; |
2251 | 2308 |
2252 dsign = bc->dsign; | 2309 dsign = bc->dsign; |
2253 nd = bc->nd; | 2310 nd = bc->nd; |
2254 nd0 = bc->nd0; | 2311 nd0 = bc->nd0; |
2255 p5 = nd + bc->e0 - 1; | 2312 p5 = nd + bc->e0 - 1; |
2256 dd = speccase = 0; | 2313 dd = speccase = 0; |
(...skipping 111 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2368 if (dd) | 2425 if (dd) |
2369 goto ret; | 2426 goto ret; |
2370 if (!b->x[0] && b->wds == 1) { | 2427 if (!b->x[0] && b->wds == 1) { |
2371 if (i < nd) | 2428 if (i < nd) |
2372 dd = 1; | 2429 dd = 1; |
2373 goto ret; | 2430 goto ret; |
2374 } | 2431 } |
2375 b = multadd(b, 10, 0); | 2432 b = multadd(b, 10, 0); |
2376 dig = quorem(b,d); | 2433 dig = quorem(b,d); |
2377 } | 2434 } |
2378 » if (b->x[0] || b->wds > 1) | 2435 » if (dig > 0 || b->x[0] || b->wds > 1) |
2379 dd = -1; | 2436 dd = -1; |
2380 ret: | 2437 ret: |
2381 Bfree(b); | 2438 Bfree(b); |
2382 Bfree(d); | 2439 Bfree(d); |
2383 #ifdef Honor_FLT_ROUNDS | 2440 #ifdef Honor_FLT_ROUNDS |
2384 if (bc->rounding != 1) { | 2441 if (bc->rounding != 1) { |
2385 if (dd < 0) { | 2442 if (dd < 0) { |
2386 if (bc->rounding == 0) { | 2443 if (bc->rounding == 0) { |
2387 if (!dsign) | 2444 if (!dsign) |
2388 goto retlow1; | 2445 goto retlow1; |
2389 } | 2446 } |
2390 else if (dsign) | 2447 else if (dsign) |
2391 goto rethi1; | 2448 goto rethi1; |
2392 } | 2449 } |
2393 else if (dd > 0) { | 2450 else if (dd > 0) { |
2394 if (bc->rounding == 0) { | 2451 if (bc->rounding == 0) { |
2395 if (dsign) | 2452 if (dsign) |
2396 goto rethi1; | 2453 goto rethi1; |
2397 goto ret1; | 2454 goto ret1; |
2398 } | 2455 } |
2399 if (!dsign) | 2456 if (!dsign) |
2400 goto rethi1; | 2457 goto rethi1; |
2401 » » » dval(rv) += 2.*ulp(rv); | 2458 » » » dval(rv) += 2.*sulp(rv,bc); |
2402 } | 2459 } |
2403 else { | 2460 else { |
2404 bc->inexact = 0; | 2461 bc->inexact = 0; |
2405 if (dsign) | 2462 if (dsign) |
2406 goto rethi1; | 2463 goto rethi1; |
2407 } | 2464 } |
2408 } | 2465 } |
2409 else | 2466 else |
2410 #endif | 2467 #endif |
2411 if (speccase) { | 2468 if (speccase) { |
2412 if (dd <= 0) | 2469 if (dd <= 0) |
2413 rv->d = 0.; | 2470 rv->d = 0.; |
2414 } | 2471 } |
2415 else if (dd < 0) { | 2472 else if (dd < 0) { |
2416 if (!dsign) /* does not happen for round-near */ | 2473 if (!dsign) /* does not happen for round-near */ |
2417 retlow1: | 2474 retlow1: |
2418 » » » dval(rv) -= ulp(rv); | 2475 » » » dval(rv) -= sulp(rv,bc); |
2419 } | 2476 } |
2420 else if (dd > 0) { | 2477 else if (dd > 0) { |
2421 if (dsign) { | 2478 if (dsign) { |
2422 rethi1: | 2479 rethi1: |
2423 » » » dval(rv) += ulp(rv); | 2480 » » » dval(rv) += sulp(rv,bc); |
2424 } | 2481 } |
2425 } | 2482 } |
2426 else { | 2483 else { |
2427 /* Exact half-way case: apply round-even rule. */ | 2484 /* Exact half-way case: apply round-even rule. */ |
2428 » » if (word1(rv) & 1) { | 2485 » » if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0
) { |
| 2486 » » » i = 1 - j; |
| 2487 » » » if (i <= 31) { |
| 2488 » » » » if (word1(rv) & (0x1 << i)) |
| 2489 » » » » » goto odd; |
| 2490 » » » » } |
| 2491 » » » else if (word0(rv) & (0x1 << (i-32))) |
| 2492 » » » » goto odd; |
| 2493 » » » } |
| 2494 » » else if (word1(rv) & 1) { |
| 2495 odd: |
2429 if (dsign) | 2496 if (dsign) |
2430 goto rethi1; | 2497 goto rethi1; |
2431 goto retlow1; | 2498 goto retlow1; |
2432 } | 2499 } |
2433 } | 2500 } |
2434 | 2501 |
2435 #ifdef Honor_FLT_ROUNDS | 2502 #ifdef Honor_FLT_ROUNDS |
2436 ret1: | 2503 ret1: |
2437 #endif | 2504 #endif |
2438 return; | 2505 return; |
2439 } | 2506 } |
2440 #endif /* NO_STRTOD_BIGCOMP */ | 2507 #endif /* NO_STRTOD_BIGCOMP */ |
2441 | 2508 |
2442 double | 2509 double |
2443 strtod | 2510 strtod |
2444 #ifdef KR_headers | 2511 #ifdef KR_headers |
2445 (s00, se) CONST char *s00; char **se; | 2512 (s00, se) CONST char *s00; char **se; |
2446 #else | 2513 #else |
2447 » (CONST char *s00, char **se) | 2514 » (const char *s00, char **se) |
2448 #endif | 2515 #endif |
2449 { | 2516 { |
2450 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1; | 2517 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1; |
2451 » int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; | 2518 » int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign; |
2452 CONST char *s, *s0, *s1; | 2519 CONST char *s, *s0, *s1; |
2453 double aadj, aadj1; | 2520 double aadj, aadj1; |
2454 Long L; | 2521 Long L; |
2455 U aadj2, adj, rv, rv0; | 2522 U aadj2, adj, rv, rv0; |
2456 ULong y, z; | 2523 ULong y, z; |
2457 BCinfo bc; | 2524 BCinfo bc; |
2458 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; | 2525 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; |
| 2526 #ifdef Avoid_Underflow |
| 2527 ULong Lsb, Lsb1; |
| 2528 #endif |
2459 #ifdef SET_INEXACT | 2529 #ifdef SET_INEXACT |
2460 int oldinexact; | 2530 int oldinexact; |
2461 #endif | 2531 #endif |
| 2532 #ifndef NO_STRTOD_BIGCOMP |
| 2533 int req_bigcomp = 0; |
| 2534 #endif |
2462 #ifdef Honor_FLT_ROUNDS /*{*/ | 2535 #ifdef Honor_FLT_ROUNDS /*{*/ |
2463 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ | 2536 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ |
2464 bc.rounding = Flt_Rounds; | 2537 bc.rounding = Flt_Rounds; |
2465 #else /*}{*/ | 2538 #else /*}{*/ |
2466 bc.rounding = 1; | 2539 bc.rounding = 1; |
2467 switch(fegetround()) { | 2540 switch(fegetround()) { |
2468 case FE_TOWARDZERO: bc.rounding = 0; break; | 2541 case FE_TOWARDZERO: bc.rounding = 0; break; |
2469 case FE_UPWARD: bc.rounding = 2; break; | 2542 case FE_UPWARD: bc.rounding = 2; break; |
2470 case FE_DOWNWARD: bc.rounding = 3; | 2543 case FE_DOWNWARD: bc.rounding = 3; |
2471 } | 2544 } |
2472 #endif /*}}*/ | 2545 #endif /*}}*/ |
2473 #endif /*}*/ | 2546 #endif /*}*/ |
2474 #ifdef USE_LOCALE | 2547 #ifdef USE_LOCALE |
2475 CONST char *s2; | 2548 CONST char *s2; |
2476 #endif | 2549 #endif |
2477 | 2550 |
2478 » sign = nz0 = nz = bc.dplen = bc.uflchk = 0; | 2551 » sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0; |
2479 dval(&rv) = 0.; | 2552 dval(&rv) = 0.; |
2480 for(s = s00;;s++) switch(*s) { | 2553 for(s = s00;;s++) switch(*s) { |
2481 case '-': | 2554 case '-': |
2482 sign = 1; | 2555 sign = 1; |
2483 /* no break */ | 2556 /* no break */ |
2484 case '+': | 2557 case '+': |
2485 if (*++s) | 2558 if (*++s) |
2486 goto break2; | 2559 goto break2; |
2487 /* no break */ | 2560 /* no break */ |
2488 case 0: | 2561 case 0: |
(...skipping 25 matching lines...) Expand all Loading... |
2514 nz0 = 1; | 2587 nz0 = 1; |
2515 while(*++s == '0') ; | 2588 while(*++s == '0') ; |
2516 if (!*s) | 2589 if (!*s) |
2517 goto ret; | 2590 goto ret; |
2518 } | 2591 } |
2519 s0 = s; | 2592 s0 = s; |
2520 y = z = 0; | 2593 y = z = 0; |
2521 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) | 2594 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) |
2522 if (nd < 9) | 2595 if (nd < 9) |
2523 y = 10*y + c - '0'; | 2596 y = 10*y + c - '0'; |
2524 » » else if (nd < 16) | 2597 » » else if (nd < DBL_DIG + 2) |
2525 z = 10*z + c - '0'; | 2598 z = 10*z + c - '0'; |
2526 nd0 = nd; | 2599 nd0 = nd; |
2527 bc.dp0 = bc.dp1 = s - s0; | 2600 bc.dp0 = bc.dp1 = s - s0; |
| 2601 for(s1 = s; s1 > s0 && *--s1 == '0'; ) |
| 2602 ++nz1; |
2528 #ifdef USE_LOCALE | 2603 #ifdef USE_LOCALE |
2529 s1 = localeconv()->decimal_point; | 2604 s1 = localeconv()->decimal_point; |
2530 if (c == *s1) { | 2605 if (c == *s1) { |
2531 c = '.'; | 2606 c = '.'; |
2532 if (*++s1) { | 2607 if (*++s1) { |
2533 s2 = s; | 2608 s2 = s; |
2534 for(;;) { | 2609 for(;;) { |
2535 if (*++s2 != *s1) { | 2610 if (*++s2 != *s1) { |
2536 c = 0; | 2611 c = 0; |
2537 break; | 2612 break; |
2538 } | 2613 } |
2539 if (!*++s1) { | 2614 if (!*++s1) { |
2540 s = s2; | 2615 s = s2; |
2541 break; | 2616 break; |
2542 } | 2617 } |
2543 } | 2618 } |
2544 } | 2619 } |
2545 } | 2620 } |
2546 #endif | 2621 #endif |
2547 if (c == '.') { | 2622 if (c == '.') { |
2548 c = *++s; | 2623 c = *++s; |
2549 bc.dp1 = s - s0; | 2624 bc.dp1 = s - s0; |
2550 bc.dplen = bc.dp1 - bc.dp0; | 2625 bc.dplen = bc.dp1 - bc.dp0; |
2551 if (!nd) { | 2626 if (!nd) { |
2552 for(; c == '0'; c = *++s) | 2627 for(; c == '0'; c = *++s) |
2553 nz++; | 2628 nz++; |
2554 if (c > '0' && c <= '9') { | 2629 if (c > '0' && c <= '9') { |
| 2630 bc.dp0 = s0 - s; |
| 2631 bc.dp1 = bc.dp0 + bc.dplen; |
2555 s0 = s; | 2632 s0 = s; |
2556 nf += nz; | 2633 nf += nz; |
2557 nz = 0; | 2634 nz = 0; |
2558 goto have_dig; | 2635 goto have_dig; |
2559 } | 2636 } |
2560 goto dig_done; | 2637 goto dig_done; |
2561 } | 2638 } |
2562 for(; c >= '0' && c <= '9'; c = *++s) { | 2639 for(; c >= '0' && c <= '9'; c = *++s) { |
2563 have_dig: | 2640 have_dig: |
2564 nz++; | 2641 nz++; |
2565 if (c -= '0') { | 2642 if (c -= '0') { |
2566 nf += nz; | 2643 nf += nz; |
2567 for(i = 1; i < nz; i++) | 2644 for(i = 1; i < nz; i++) |
2568 if (nd++ < 9) | 2645 if (nd++ < 9) |
2569 y *= 10; | 2646 y *= 10; |
2570 » » » » » else if (nd <= DBL_DIG + 1) | 2647 » » » » » else if (nd <= DBL_DIG + 2) |
2571 z *= 10; | 2648 z *= 10; |
2572 if (nd++ < 9) | 2649 if (nd++ < 9) |
2573 y = 10*y + c; | 2650 y = 10*y + c; |
2574 » » » » else if (nd <= DBL_DIG + 1) | 2651 » » » » else if (nd <= DBL_DIG + 2) |
2575 z = 10*z + c; | 2652 z = 10*z + c; |
2576 » » » » nz = 0; | 2653 » » » » nz = nz1 = 0; |
2577 } | 2654 } |
2578 } | 2655 } |
2579 } | 2656 } |
2580 dig_done: | 2657 dig_done: |
2581 e = 0; | 2658 e = 0; |
2582 if (c == 'e' || c == 'E') { | 2659 if (c == 'e' || c == 'E') { |
2583 if (!nd && !nz && !nz0) { | 2660 if (!nd && !nz && !nz0) { |
2584 goto ret0; | 2661 goto ret0; |
2585 } | 2662 } |
2586 s00 = s; | 2663 s00 = s; |
(...skipping 69 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2656 } | 2733 } |
2657 bc.e0 = e1 = e -= nf; | 2734 bc.e0 = e1 = e -= nf; |
2658 | 2735 |
2659 /* Now we have nd0 digits, starting at s0, followed by a | 2736 /* Now we have nd0 digits, starting at s0, followed by a |
2660 * decimal point, followed by nd-nd0 digits. The number we're | 2737 * decimal point, followed by nd-nd0 digits. The number we're |
2661 * after is the integer represented by those digits times | 2738 * after is the integer represented by those digits times |
2662 * 10**e */ | 2739 * 10**e */ |
2663 | 2740 |
2664 if (!nd0) | 2741 if (!nd0) |
2665 nd0 = nd; | 2742 nd0 = nd; |
2666 » k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; | 2743 » k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2; |
2667 dval(&rv) = y; | 2744 dval(&rv) = y; |
2668 if (k > 9) { | 2745 if (k > 9) { |
2669 #ifdef SET_INEXACT | 2746 #ifdef SET_INEXACT |
2670 if (k > DBL_DIG) | 2747 if (k > DBL_DIG) |
2671 oldinexact = get_inexact(); | 2748 oldinexact = get_inexact(); |
2672 #endif | 2749 #endif |
2673 dval(&rv) = tens[k - 9] * dval(&rv) + z; | 2750 dval(&rv) = tens[k - 9] * dval(&rv) + z; |
2674 } | 2751 } |
2675 bd0 = 0; | 2752 bd0 = 0; |
2676 if (nd <= DBL_DIG | 2753 if (nd <= DBL_DIG |
2677 #ifndef RND_PRODQUOT | 2754 #ifndef RND_PRODQUOT |
2678 #ifndef Honor_FLT_ROUNDS | 2755 #ifndef Honor_FLT_ROUNDS |
2679 && Flt_Rounds == 1 | 2756 && Flt_Rounds == 1 |
2680 #endif | 2757 #endif |
2681 #endif | 2758 #endif |
2682 ) { | 2759 ) { |
2683 if (!e) | 2760 if (!e) |
2684 goto ret; | 2761 goto ret; |
| 2762 #ifndef ROUND_BIASED_without_Round_Up |
2685 if (e > 0) { | 2763 if (e > 0) { |
2686 if (e <= Ten_pmax) { | 2764 if (e <= Ten_pmax) { |
2687 #ifdef VAX | 2765 #ifdef VAX |
2688 goto vax_ovfl_check; | 2766 goto vax_ovfl_check; |
2689 #else | 2767 #else |
2690 #ifdef Honor_FLT_ROUNDS | 2768 #ifdef Honor_FLT_ROUNDS |
2691 /* round correctly FLT_ROUNDS = 2 or 3 */ | 2769 /* round correctly FLT_ROUNDS = 2 or 3 */ |
2692 if (sign) { | 2770 if (sign) { |
2693 rv.d = -rv.d; | 2771 rv.d = -rv.d; |
2694 sign = 0; | 2772 sign = 0; |
(...skipping 40 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2735 /* round correctly FLT_ROUNDS = 2 or 3 */ | 2813 /* round correctly FLT_ROUNDS = 2 or 3 */ |
2736 if (sign) { | 2814 if (sign) { |
2737 rv.d = -rv.d; | 2815 rv.d = -rv.d; |
2738 sign = 0; | 2816 sign = 0; |
2739 } | 2817 } |
2740 #endif | 2818 #endif |
2741 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); | 2819 /* rv = */ rounded_quotient(dval(&rv), tens[-e]); |
2742 goto ret; | 2820 goto ret; |
2743 } | 2821 } |
2744 #endif | 2822 #endif |
| 2823 #endif /* ROUND_BIASED_without_Round_Up */ |
2745 } | 2824 } |
2746 e1 += nd - k; | 2825 e1 += nd - k; |
2747 | 2826 |
2748 #ifdef IEEE_Arith | 2827 #ifdef IEEE_Arith |
2749 #ifdef SET_INEXACT | 2828 #ifdef SET_INEXACT |
2750 bc.inexact = 1; | 2829 bc.inexact = 1; |
2751 if (k <= DBL_DIG) | 2830 if (k <= DBL_DIG) |
2752 oldinexact = get_inexact(); | 2831 oldinexact = get_inexact(); |
2753 #endif | 2832 #endif |
2754 #ifdef Avoid_Underflow | 2833 #ifdef Avoid_Underflow |
(...skipping 12 matching lines...) Expand all Loading... |
2767 | 2846 |
2768 /* Get starting approximation = rv * 10**e1 */ | 2847 /* Get starting approximation = rv * 10**e1 */ |
2769 | 2848 |
2770 if (e1 > 0) { | 2849 if (e1 > 0) { |
2771 i = e1 & 15; | 2850 i = e1 & 15; |
2772 if (i) | 2851 if (i) |
2773 dval(&rv) *= tens[i]; | 2852 dval(&rv) *= tens[i]; |
2774 if (e1 &= ~15) { | 2853 if (e1 &= ~15) { |
2775 if (e1 > DBL_MAX_10_EXP) { | 2854 if (e1 > DBL_MAX_10_EXP) { |
2776 ovfl: | 2855 ovfl: |
2777 #ifndef NO_ERRNO | |
2778 errno = ERANGE; | |
2779 #endif | |
2780 /* Can't trust HUGE_VAL */ | 2856 /* Can't trust HUGE_VAL */ |
2781 #ifdef IEEE_Arith | 2857 #ifdef IEEE_Arith |
2782 #ifdef Honor_FLT_ROUNDS | 2858 #ifdef Honor_FLT_ROUNDS |
2783 switch(bc.rounding) { | 2859 switch(bc.rounding) { |
2784 case 0: /* toward 0 */ | 2860 case 0: /* toward 0 */ |
2785 case 3: /* toward -infinity */ | 2861 case 3: /* toward -infinity */ |
2786 word0(&rv) = Big0; | 2862 word0(&rv) = Big0; |
2787 word1(&rv) = Big1; | 2863 word1(&rv) = Big1; |
2788 break; | 2864 break; |
2789 default: | 2865 default: |
2790 word0(&rv) = Exp_mask; | 2866 word0(&rv) = Exp_mask; |
2791 word1(&rv) = 0; | 2867 word1(&rv) = 0; |
2792 } | 2868 } |
2793 #else /*Honor_FLT_ROUNDS*/ | 2869 #else /*Honor_FLT_ROUNDS*/ |
2794 word0(&rv) = Exp_mask; | 2870 word0(&rv) = Exp_mask; |
2795 word1(&rv) = 0; | 2871 word1(&rv) = 0; |
2796 #endif /*Honor_FLT_ROUNDS*/ | 2872 #endif /*Honor_FLT_ROUNDS*/ |
2797 #ifdef SET_INEXACT | 2873 #ifdef SET_INEXACT |
2798 /* set overflow bit */ | 2874 /* set overflow bit */ |
2799 dval(&rv0) = 1e300; | 2875 dval(&rv0) = 1e300; |
2800 dval(&rv0) *= dval(&rv0); | 2876 dval(&rv0) *= dval(&rv0); |
2801 #endif | 2877 #endif |
2802 #else /*IEEE_Arith*/ | 2878 #else /*IEEE_Arith*/ |
2803 word0(&rv) = Big0; | 2879 word0(&rv) = Big0; |
2804 word1(&rv) = Big1; | 2880 word1(&rv) = Big1; |
2805 #endif /*IEEE_Arith*/ | 2881 #endif /*IEEE_Arith*/ |
| 2882 range_err: |
| 2883 if (bd0) { |
| 2884 Bfree(bb); |
| 2885 Bfree(bd); |
| 2886 Bfree(bs); |
| 2887 Bfree(bd0); |
| 2888 Bfree(delta); |
| 2889 } |
| 2890 #ifndef NO_ERRNO |
| 2891 errno = ERANGE; |
| 2892 #endif |
2806 goto ret; | 2893 goto ret; |
2807 } | 2894 } |
2808 e1 >>= 4; | 2895 e1 >>= 4; |
2809 for(j = 0; e1 > 1; j++, e1 >>= 1) | 2896 for(j = 0; e1 > 1; j++, e1 >>= 1) |
2810 if (e1 & 1) | 2897 if (e1 & 1) |
2811 dval(&rv) *= bigtens[j]; | 2898 dval(&rv) *= bigtens[j]; |
2812 /* The last multiplication could overflow. */ | 2899 /* The last multiplication could overflow. */ |
2813 word0(&rv) -= P*Exp_msk1; | 2900 word0(&rv) -= P*Exp_msk1; |
2814 dval(&rv) *= bigtens[j]; | 2901 dval(&rv) *= bigtens[j]; |
2815 if ((z = word0(&rv) & Exp_mask) | 2902 if ((z = word0(&rv) & Exp_mask) |
(...skipping 20 matching lines...) Expand all Loading... |
2836 #ifdef Avoid_Underflow | 2923 #ifdef Avoid_Underflow |
2837 if (e1 & Scale_Bit) | 2924 if (e1 & Scale_Bit) |
2838 bc.scale = 2*P; | 2925 bc.scale = 2*P; |
2839 for(j = 0; e1 > 0; j++, e1 >>= 1) | 2926 for(j = 0; e1 > 0; j++, e1 >>= 1) |
2840 if (e1 & 1) | 2927 if (e1 & 1) |
2841 dval(&rv) *= tinytens[j]; | 2928 dval(&rv) *= tinytens[j]; |
2842 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) | 2929 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) |
2843 >> Exp_shift)) > 0) { | 2930 >> Exp_shift)) > 0) { |
2844 /* scaled rv is denormal; clear j low bits */ | 2931 /* scaled rv is denormal; clear j low bits */ |
2845 if (j >= 32) { | 2932 if (j >= 32) { |
| 2933 if (j > 54) |
| 2934 goto undfl; |
2846 word1(&rv) = 0; | 2935 word1(&rv) = 0; |
2847 if (j >= 53) | 2936 if (j >= 53) |
2848 word0(&rv) = (P+2)*Exp_msk1; | 2937 word0(&rv) = (P+2)*Exp_msk1; |
2849 else | 2938 else |
2850 word0(&rv) &= 0xffffffff << (j-32); | 2939 word0(&rv) &= 0xffffffff << (j-32); |
2851 } | 2940 } |
2852 else | 2941 else |
2853 word1(&rv) &= 0xffffffff << j; | 2942 word1(&rv) &= 0xffffffff << j; |
2854 } | 2943 } |
2855 #else | 2944 #else |
2856 for(j = 0; e1 > 1; j++, e1 >>= 1) | 2945 for(j = 0; e1 > 1; j++, e1 >>= 1) |
2857 if (e1 & 1) | 2946 if (e1 & 1) |
2858 dval(&rv) *= tinytens[j]; | 2947 dval(&rv) *= tinytens[j]; |
2859 /* The last multiplication could underflow. */ | 2948 /* The last multiplication could underflow. */ |
2860 dval(&rv0) = dval(&rv); | 2949 dval(&rv0) = dval(&rv); |
2861 dval(&rv) *= tinytens[j]; | 2950 dval(&rv) *= tinytens[j]; |
2862 if (!dval(&rv)) { | 2951 if (!dval(&rv)) { |
2863 dval(&rv) = 2.*dval(&rv0); | 2952 dval(&rv) = 2.*dval(&rv0); |
2864 dval(&rv) *= tinytens[j]; | 2953 dval(&rv) *= tinytens[j]; |
2865 #endif | 2954 #endif |
2866 if (!dval(&rv)) { | 2955 if (!dval(&rv)) { |
2867 undfl: | 2956 undfl: |
2868 dval(&rv) = 0.; | 2957 dval(&rv) = 0.; |
2869 #ifndef NO_ERRNO | 2958 » » » » » goto range_err; |
2870 » » » » » errno = ERANGE; | |
2871 #endif | |
2872 » » » » » goto ret; | |
2873 } | 2959 } |
2874 #ifndef Avoid_Underflow | 2960 #ifndef Avoid_Underflow |
2875 word0(&rv) = Tiny0; | 2961 word0(&rv) = Tiny0; |
2876 word1(&rv) = Tiny1; | 2962 word1(&rv) = Tiny1; |
2877 /* The refinement below will clean | 2963 /* The refinement below will clean |
2878 * this approximation up. | 2964 * this approximation up. |
2879 */ | 2965 */ |
2880 } | 2966 } |
2881 #endif | 2967 #endif |
2882 } | 2968 } |
2883 } | 2969 } |
2884 | 2970 |
2885 /* Now the hard part -- adjusting rv to the correct value.*/ | 2971 /* Now the hard part -- adjusting rv to the correct value.*/ |
2886 | 2972 |
2887 /* Put digits into bd: true value = bd * 10^e */ | 2973 /* Put digits into bd: true value = bd * 10^e */ |
2888 | 2974 |
2889 » bc.nd = nd; | 2975 » bc.nd = nd - nz1; |
2890 #ifndef NO_STRTOD_BIGCOMP | 2976 #ifndef NO_STRTOD_BIGCOMP |
2891 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */ | 2977 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */ |
2892 /* to silence an erroneous warning about bc.nd0 */ | 2978 /* to silence an erroneous warning about bc.nd0 */ |
2893 /* possibly not being initialized. */ | 2979 /* possibly not being initialized. */ |
2894 if (nd > strtod_diglim) { | 2980 if (nd > strtod_diglim) { |
2895 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */ | 2981 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */ |
2896 /* minimum number of decimal digits to distinguish double values
*/ | 2982 /* minimum number of decimal digits to distinguish double values
*/ |
2897 /* in IEEE arithmetic. */ | 2983 /* in IEEE arithmetic. */ |
2898 i = j = 18; | 2984 i = j = 18; |
2899 if (i > nd0) | 2985 if (i > nd0) |
2900 j += bc.dplen; | 2986 j += bc.dplen; |
2901 for(;;) { | 2987 for(;;) { |
2902 » » » if (--j <= bc.dp1 && j >= bc.dp0) | 2988 » » » if (--j < bc.dp1 && j >= bc.dp0) |
2903 j = bc.dp0 - 1; | 2989 j = bc.dp0 - 1; |
2904 if (s0[j] != '0') | 2990 if (s0[j] != '0') |
2905 break; | 2991 break; |
2906 --i; | 2992 --i; |
2907 } | 2993 } |
2908 e += nd - i; | 2994 e += nd - i; |
2909 nd = i; | 2995 nd = i; |
2910 if (nd0 > nd) | 2996 if (nd0 > nd) |
2911 nd0 = nd; | 2997 nd0 = nd; |
2912 if (nd < 9) { /* must recompute y */ | 2998 if (nd < 9) { /* must recompute y */ |
(...skipping 24 matching lines...) Expand all Loading... |
2937 if (bbe >= 0) | 3023 if (bbe >= 0) |
2938 bb2 += bbe; | 3024 bb2 += bbe; |
2939 else | 3025 else |
2940 bd2 -= bbe; | 3026 bd2 -= bbe; |
2941 bs2 = bb2; | 3027 bs2 = bb2; |
2942 #ifdef Honor_FLT_ROUNDS | 3028 #ifdef Honor_FLT_ROUNDS |
2943 if (bc.rounding != 1) | 3029 if (bc.rounding != 1) |
2944 bs2++; | 3030 bs2++; |
2945 #endif | 3031 #endif |
2946 #ifdef Avoid_Underflow | 3032 #ifdef Avoid_Underflow |
| 3033 Lsb = LSB; |
| 3034 Lsb1 = 0; |
2947 j = bbe - bc.scale; | 3035 j = bbe - bc.scale; |
2948 i = j + bbbits - 1; /* logb(rv) */ | 3036 i = j + bbbits - 1; /* logb(rv) */ |
2949 » » if (i < Emin)» /* denormal */ | 3037 » » j = P + 1 - bbbits; |
2950 » » » j += P - Emin; | 3038 » » if (i < Emin) {»/* denormal */ |
2951 » » else | 3039 » » » i = Emin - i; |
2952 » » » j = P + 1 - bbbits; | 3040 » » » j -= i; |
| 3041 » » » if (i < 32) |
| 3042 » » » » Lsb <<= i; |
| 3043 » » » else if (i < 52) |
| 3044 » » » » Lsb1 = Lsb << (i-32); |
| 3045 » » » else |
| 3046 » » » » Lsb1 = Exp_mask; |
| 3047 » » » } |
2953 #else /*Avoid_Underflow*/ | 3048 #else /*Avoid_Underflow*/ |
2954 #ifdef Sudden_Underflow | 3049 #ifdef Sudden_Underflow |
2955 #ifdef IBM | 3050 #ifdef IBM |
2956 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); | 3051 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); |
2957 #else | 3052 #else |
2958 j = P + 1 - bbbits; | 3053 j = P + 1 - bbbits; |
2959 #endif | 3054 #endif |
2960 #else /*Sudden_Underflow*/ | 3055 #else /*Sudden_Underflow*/ |
2961 j = bbe; | 3056 j = bbe; |
2962 i = j + bbbits - 1; /* logb(rv) */ | 3057 i = j + bbbits - 1; /* logb(rv) */ |
(...skipping 27 matching lines...) Expand all Loading... |
2990 if (bd5 > 0) | 3085 if (bd5 > 0) |
2991 bd = pow5mult(bd, bd5); | 3086 bd = pow5mult(bd, bd5); |
2992 if (bd2 > 0) | 3087 if (bd2 > 0) |
2993 bd = lshift(bd, bd2); | 3088 bd = lshift(bd, bd2); |
2994 if (bs2 > 0) | 3089 if (bs2 > 0) |
2995 bs = lshift(bs, bs2); | 3090 bs = lshift(bs, bs2); |
2996 delta = diff(bb, bd); | 3091 delta = diff(bb, bd); |
2997 bc.dsign = delta->sign; | 3092 bc.dsign = delta->sign; |
2998 delta->sign = 0; | 3093 delta->sign = 0; |
2999 i = cmp(delta, bs); | 3094 i = cmp(delta, bs); |
3000 #ifndef NO_STRTOD_BIGCOMP | 3095 #ifndef NO_STRTOD_BIGCOMP /*{*/ |
3001 if (bc.nd > nd && i <= 0) { | 3096 if (bc.nd > nd && i <= 0) { |
3002 » » » if (bc.dsign) | 3097 » » » if (bc.dsign) { |
3003 » » » » break;» /* Must use bigcomp(). */ | 3098 » » » » /* Must use bigcomp(). */ |
| 3099 » » » » req_bigcomp = 1; |
| 3100 » » » » break; |
| 3101 » » » » } |
3004 #ifdef Honor_FLT_ROUNDS | 3102 #ifdef Honor_FLT_ROUNDS |
3005 if (bc.rounding != 1) { | 3103 if (bc.rounding != 1) { |
3006 » » » » if (i < 0) | 3104 » » » » if (i < 0) { |
| 3105 » » » » » req_bigcomp = 1; |
3007 break; | 3106 break; |
| 3107 } |
3008 } | 3108 } |
3009 else | 3109 else |
3010 #endif | 3110 #endif |
3011 { | |
3012 bc.nd = nd; | |
3013 i = -1; /* Discarded digits make delta smaller.
*/ | 3111 i = -1; /* Discarded digits make delta smaller.
*/ |
3014 } | |
3015 } | 3112 } |
3016 #endif | 3113 #endif /*}*/ |
3017 #ifdef Honor_FLT_ROUNDS | 3114 #ifdef Honor_FLT_ROUNDS /*{*/ |
3018 if (bc.rounding != 1) { | 3115 if (bc.rounding != 1) { |
3019 if (i < 0) { | 3116 if (i < 0) { |
3020 /* Error is less than an ulp */ | 3117 /* Error is less than an ulp */ |
3021 if (!delta->x[0] && delta->wds <= 1) { | 3118 if (!delta->x[0] && delta->wds <= 1) { |
3022 /* exact */ | 3119 /* exact */ |
3023 #ifdef SET_INEXACT | 3120 #ifdef SET_INEXACT |
3024 bc.inexact = 0; | 3121 bc.inexact = 0; |
3025 #endif | 3122 #endif |
3026 break; | 3123 break; |
3027 } | 3124 } |
(...skipping 13 matching lines...) Expand all Loading... |
3041 #else | 3138 #else |
3042 if (y) | 3139 if (y) |
3043 #endif | 3140 #endif |
3044 { | 3141 { |
3045 delta = lshift(delta,Log2P); | 3142 delta = lshift(delta,Log2P); |
3046 if (cmp(delta, bs) <= 0) | 3143 if (cmp(delta, bs) <= 0) |
3047 adj.d = -0.5; | 3144 adj.d = -0.5; |
3048 } | 3145 } |
3049 } | 3146 } |
3050 apply_adj: | 3147 apply_adj: |
3051 #ifdef Avoid_Underflow | 3148 #ifdef Avoid_Underflow /*{*/ |
3052 if (bc.scale && (y = word0(&rv) & Exp_ma
sk) | 3149 if (bc.scale && (y = word0(&rv) & Exp_ma
sk) |
3053 <= 2*P*Exp_msk1) | 3150 <= 2*P*Exp_msk1) |
3054 word0(&adj) += (2*P+1)*Exp_msk1 - y; | 3151 word0(&adj) += (2*P+1)*Exp_msk1 - y; |
3055 #else | 3152 #else |
3056 #ifdef Sudden_Underflow | 3153 #ifdef Sudden_Underflow |
3057 if ((word0(&rv) & Exp_mask) <= | 3154 if ((word0(&rv) & Exp_mask) <= |
3058 P*Exp_msk1) { | 3155 P*Exp_msk1) { |
3059 word0(&rv) += P*Exp_msk1; | 3156 word0(&rv) += P*Exp_msk1; |
3060 dval(&rv) += adj.d*ulp(dval(&rv)
); | 3157 dval(&rv) += adj.d*ulp(dval(&rv)
); |
3061 word0(&rv) -= P*Exp_msk1; | 3158 word0(&rv) -= P*Exp_msk1; |
3062 } | 3159 } |
3063 else | 3160 else |
3064 #endif /*Sudden_Underflow*/ | 3161 #endif /*Sudden_Underflow*/ |
3065 #endif /*Avoid_Underflow*/ | 3162 #endif /*Avoid_Underflow}*/ |
3066 dval(&rv) += adj.d*ulp(&rv); | 3163 dval(&rv) += adj.d*ulp(&rv); |
3067 } | 3164 } |
3068 break; | 3165 break; |
3069 } | 3166 } |
3070 adj.d = ratio(delta, bs); | 3167 adj.d = ratio(delta, bs); |
3071 if (adj.d < 1.) | 3168 if (adj.d < 1.) |
3072 adj.d = 1.; | 3169 adj.d = 1.; |
3073 if (adj.d <= 0x7ffffffe) { | 3170 if (adj.d <= 0x7ffffffe) { |
3074 /* adj = rounding ? ceil(adj) : floor(adj); */ | 3171 /* adj = rounding ? ceil(adj) : floor(adj); */ |
3075 y = adj.d; | 3172 y = adj.d; |
3076 if (y != adj.d) { | 3173 if (y != adj.d) { |
3077 if (!((bc.rounding>>1) ^ bc.dsign)) | 3174 if (!((bc.rounding>>1) ^ bc.dsign)) |
3078 y++; | 3175 y++; |
3079 adj.d = y; | 3176 adj.d = y; |
3080 } | 3177 } |
3081 } | 3178 } |
3082 #ifdef Avoid_Underflow | 3179 #ifdef Avoid_Underflow /*{*/ |
3083 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_m
sk1) | 3180 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_m
sk1) |
3084 word0(&adj) += (2*P+1)*Exp_msk1 - y; | 3181 word0(&adj) += (2*P+1)*Exp_msk1 - y; |
3085 #else | 3182 #else |
3086 #ifdef Sudden_Underflow | 3183 #ifdef Sudden_Underflow |
3087 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { | 3184 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { |
3088 word0(&rv) += P*Exp_msk1; | 3185 word0(&rv) += P*Exp_msk1; |
3089 adj.d *= ulp(dval(&rv)); | 3186 adj.d *= ulp(dval(&rv)); |
3090 if (bc.dsign) | 3187 if (bc.dsign) |
3091 dval(&rv) += adj.d; | 3188 dval(&rv) += adj.d; |
3092 else | 3189 else |
3093 dval(&rv) -= adj.d; | 3190 dval(&rv) -= adj.d; |
3094 word0(&rv) -= P*Exp_msk1; | 3191 word0(&rv) -= P*Exp_msk1; |
3095 goto cont; | 3192 goto cont; |
3096 } | 3193 } |
3097 #endif /*Sudden_Underflow*/ | 3194 #endif /*Sudden_Underflow*/ |
3098 #endif /*Avoid_Underflow*/ | 3195 #endif /*Avoid_Underflow}*/ |
3099 adj.d *= ulp(&rv); | 3196 adj.d *= ulp(&rv); |
3100 if (bc.dsign) { | 3197 if (bc.dsign) { |
3101 if (word0(&rv) == Big0 && word1(&rv) == Big1) | 3198 if (word0(&rv) == Big0 && word1(&rv) == Big1) |
3102 goto ovfl; | 3199 goto ovfl; |
3103 dval(&rv) += adj.d; | 3200 dval(&rv) += adj.d; |
3104 } | 3201 } |
3105 else | 3202 else |
3106 dval(&rv) -= adj.d; | 3203 dval(&rv) -= adj.d; |
3107 goto cont; | 3204 goto cont; |
3108 } | 3205 } |
3109 #endif /*Honor_FLT_ROUNDS*/ | 3206 #endif /*}Honor_FLT_ROUNDS*/ |
3110 | 3207 |
3111 if (i < 0) { | 3208 if (i < 0) { |
3112 /* Error is less than half an ulp -- check for | 3209 /* Error is less than half an ulp -- check for |
3113 * special case of mantissa a power of two. | 3210 * special case of mantissa a power of two. |
3114 */ | 3211 */ |
3115 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask | 3212 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask |
3116 #ifdef IEEE_Arith | 3213 #ifdef IEEE_Arith /*{*/ |
3117 #ifdef Avoid_Underflow | 3214 #ifdef Avoid_Underflow |
3118 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 | 3215 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 |
3119 #else | 3216 #else |
3120 || (word0(&rv) & Exp_mask) <= Exp_msk1 | 3217 || (word0(&rv) & Exp_mask) <= Exp_msk1 |
3121 #endif | 3218 #endif |
3122 #endif | 3219 #endif /*}*/ |
3123 ) { | 3220 ) { |
3124 #ifdef SET_INEXACT | 3221 #ifdef SET_INEXACT |
3125 if (!delta->x[0] && delta->wds <= 1) | 3222 if (!delta->x[0] && delta->wds <= 1) |
3126 bc.inexact = 0; | 3223 bc.inexact = 0; |
3127 #endif | 3224 #endif |
3128 break; | 3225 break; |
3129 } | 3226 } |
3130 if (!delta->x[0] && delta->wds <= 1) { | 3227 if (!delta->x[0] && delta->wds <= 1) { |
3131 /* exact result */ | 3228 /* exact result */ |
3132 #ifdef SET_INEXACT | 3229 #ifdef SET_INEXACT |
(...skipping 10 matching lines...) Expand all Loading... |
3143 /* exactly half-way between */ | 3240 /* exactly half-way between */ |
3144 if (bc.dsign) { | 3241 if (bc.dsign) { |
3145 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 | 3242 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 |
3146 && word1(&rv) == ( | 3243 && word1(&rv) == ( |
3147 #ifdef Avoid_Underflow | 3244 #ifdef Avoid_Underflow |
3148 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1
) | 3245 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1
) |
3149 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : | 3246 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |
3150 #endif | 3247 #endif |
3151 0xffffffff)) { | 3248 0xffffffff)) { |
3152 /*boundary case -- increment exponent*/ | 3249 /*boundary case -- increment exponent*/ |
| 3250 if (word0(&rv) == Big0 && word1(&rv) ==
Big1) |
| 3251 goto ovfl; |
3153 word0(&rv) = (word0(&rv) & Exp_mask) | 3252 word0(&rv) = (word0(&rv) & Exp_mask) |
3154 + Exp_msk1 | 3253 + Exp_msk1 |
3155 #ifdef IBM | 3254 #ifdef IBM |
3156 | Exp_msk1 >> 4 | 3255 | Exp_msk1 >> 4 |
3157 #endif | 3256 #endif |
3158 ; | 3257 ; |
3159 word1(&rv) = 0; | 3258 word1(&rv) = 0; |
3160 #ifdef Avoid_Underflow | 3259 #ifdef Avoid_Underflow |
3161 bc.dsign = 0; | 3260 bc.dsign = 0; |
3162 #endif | 3261 #endif |
(...skipping 40 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3203 } | 3302 } |
3204 } | 3303 } |
3205 #endif /*Avoid_Underflow*/ | 3304 #endif /*Avoid_Underflow*/ |
3206 L = (word0(&rv) & Exp_mask) - Exp_msk1; | 3305 L = (word0(&rv) & Exp_mask) - Exp_msk1; |
3207 #endif /*Sudden_Underflow}}*/ | 3306 #endif /*Sudden_Underflow}}*/ |
3208 word0(&rv) = L | Bndry_mask1; | 3307 word0(&rv) = L | Bndry_mask1; |
3209 word1(&rv) = 0xffffffff; | 3308 word1(&rv) = 0xffffffff; |
3210 #ifdef IBM | 3309 #ifdef IBM |
3211 goto cont; | 3310 goto cont; |
3212 #else | 3311 #else |
| 3312 #ifndef NO_STRTOD_BIGCOMP |
| 3313 if (bc.nd > nd) |
| 3314 goto cont; |
| 3315 #endif |
3213 break; | 3316 break; |
3214 #endif | 3317 #endif |
3215 } | 3318 } |
3216 #ifndef ROUND_BIASED | 3319 #ifndef ROUND_BIASED |
| 3320 #ifdef Avoid_Underflow |
| 3321 if (Lsb1) { |
| 3322 if (!(word0(&rv) & Lsb1)) |
| 3323 break; |
| 3324 } |
| 3325 else if (!(word1(&rv) & Lsb)) |
| 3326 break; |
| 3327 #else |
3217 if (!(word1(&rv) & LSB)) | 3328 if (!(word1(&rv) & LSB)) |
3218 break; | 3329 break; |
3219 #endif | 3330 #endif |
| 3331 #endif |
3220 if (bc.dsign) | 3332 if (bc.dsign) |
| 3333 #ifdef Avoid_Underflow |
| 3334 dval(&rv) += sulp(&rv, &bc); |
| 3335 #else |
3221 dval(&rv) += ulp(&rv); | 3336 dval(&rv) += ulp(&rv); |
| 3337 #endif |
3222 #ifndef ROUND_BIASED | 3338 #ifndef ROUND_BIASED |
3223 else { | 3339 else { |
| 3340 #ifdef Avoid_Underflow |
| 3341 dval(&rv) -= sulp(&rv, &bc); |
| 3342 #else |
3224 dval(&rv) -= ulp(&rv); | 3343 dval(&rv) -= ulp(&rv); |
| 3344 #endif |
3225 #ifndef Sudden_Underflow | 3345 #ifndef Sudden_Underflow |
3226 if (!dval(&rv)) { | 3346 if (!dval(&rv)) { |
3227 if (bc.nd >nd) { | 3347 if (bc.nd >nd) { |
3228 bc.uflchk = 1; | 3348 bc.uflchk = 1; |
3229 break; | 3349 break; |
3230 } | 3350 } |
3231 goto undfl; | 3351 goto undfl; |
3232 } | 3352 } |
3233 #endif | 3353 #endif |
3234 } | 3354 } |
(...skipping 72 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3307 if (bc.scale && y <= 2*P*Exp_msk1) { | 3427 if (bc.scale && y <= 2*P*Exp_msk1) { |
3308 if (aadj <= 0x7fffffff) { | 3428 if (aadj <= 0x7fffffff) { |
3309 if ((z = (ULong)aadj) <= 0) | 3429 if ((z = (ULong)aadj) <= 0) |
3310 z = 1; | 3430 z = 1; |
3311 aadj = z; | 3431 aadj = z; |
3312 aadj1 = bc.dsign ? aadj : -aadj; | 3432 aadj1 = bc.dsign ? aadj : -aadj; |
3313 } | 3433 } |
3314 dval(&aadj2) = aadj1; | 3434 dval(&aadj2) = aadj1; |
3315 word0(&aadj2) += (2*P+1)*Exp_msk1 - y; | 3435 word0(&aadj2) += (2*P+1)*Exp_msk1 - y; |
3316 aadj1 = dval(&aadj2); | 3436 aadj1 = dval(&aadj2); |
| 3437 adj.d = aadj1 * ulp(&rv); |
| 3438 dval(&rv) += adj.d; |
| 3439 if (rv.d == 0.) |
| 3440 #ifdef NO_STRTOD_BIGCOMP |
| 3441 goto undfl; |
| 3442 #else |
| 3443 { |
| 3444 req_bigcomp = 1; |
| 3445 break; |
| 3446 } |
| 3447 #endif |
3317 } | 3448 } |
3318 » » » adj.d = aadj1 * ulp(&rv); | 3449 » » » else { |
3319 » » » dval(&rv) += adj.d; | 3450 » » » » adj.d = aadj1 * ulp(&rv); |
| 3451 » » » » dval(&rv) += adj.d; |
| 3452 » » » » } |
3320 #else | 3453 #else |
3321 #ifdef Sudden_Underflow | 3454 #ifdef Sudden_Underflow |
3322 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { | 3455 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) { |
3323 dval(&rv0) = dval(&rv); | 3456 dval(&rv0) = dval(&rv); |
3324 word0(&rv) += P*Exp_msk1; | 3457 word0(&rv) += P*Exp_msk1; |
3325 adj.d = aadj1 * ulp(&rv); | 3458 adj.d = aadj1 * ulp(&rv); |
3326 dval(&rv) += adj.d; | 3459 dval(&rv) += adj.d; |
3327 #ifdef IBM | 3460 #ifdef IBM |
3328 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) | 3461 if ((word0(&rv) & Exp_mask) < P*Exp_msk1) |
3329 #else | 3462 #else |
(...skipping 62 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3392 Bfree(bd); | 3525 Bfree(bd); |
3393 Bfree(bs); | 3526 Bfree(bs); |
3394 Bfree(delta); | 3527 Bfree(delta); |
3395 } | 3528 } |
3396 Bfree(bb); | 3529 Bfree(bb); |
3397 Bfree(bd); | 3530 Bfree(bd); |
3398 Bfree(bs); | 3531 Bfree(bs); |
3399 Bfree(bd0); | 3532 Bfree(bd0); |
3400 Bfree(delta); | 3533 Bfree(delta); |
3401 #ifndef NO_STRTOD_BIGCOMP | 3534 #ifndef NO_STRTOD_BIGCOMP |
3402 » if (bc.nd > nd) | 3535 » if (req_bigcomp) { |
| 3536 » » bd0 = 0; |
| 3537 » » bc.e0 += nz1; |
3403 bigcomp(&rv, s0, &bc); | 3538 bigcomp(&rv, s0, &bc); |
| 3539 y = word0(&rv) & Exp_mask; |
| 3540 if (y == Exp_mask) |
| 3541 goto ovfl; |
| 3542 if (y == 0 && rv.d == 0.) |
| 3543 goto undfl; |
| 3544 } |
3404 #endif | 3545 #endif |
3405 #ifdef SET_INEXACT | 3546 #ifdef SET_INEXACT |
3406 if (bc.inexact) { | 3547 if (bc.inexact) { |
3407 if (!oldinexact) { | 3548 if (!oldinexact) { |
3408 word0(&rv0) = Exp_1 + (70 << Exp_shift); | 3549 word0(&rv0) = Exp_1 + (70 << Exp_shift); |
3409 word1(&rv0) = 0; | 3550 word1(&rv0) = 0; |
3410 dval(&rv0) += 1.; | 3551 dval(&rv0) += 1.; |
3411 } | 3552 } |
3412 } | 3553 } |
3413 else if (!oldinexact) | 3554 else if (!oldinexact) |
(...skipping 52 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3466 #ifndef MULTIPLE_THREADS | 3607 #ifndef MULTIPLE_THREADS |
3467 dtoa_result = | 3608 dtoa_result = |
3468 #endif | 3609 #endif |
3469 (char *)(r+1); | 3610 (char *)(r+1); |
3470 } | 3611 } |
3471 | 3612 |
3472 static char * | 3613 static char * |
3473 #ifdef KR_headers | 3614 #ifdef KR_headers |
3474 nrv_alloc(s, rve, n) char *s, **rve; int n; | 3615 nrv_alloc(s, rve, n) char *s, **rve; int n; |
3475 #else | 3616 #else |
3476 nrv_alloc(CONST char *s, char **rve, int n) | 3617 nrv_alloc(const char *s, char **rve, int n) |
3477 #endif | 3618 #endif |
3478 { | 3619 { |
3479 char *rv, *t; | 3620 char *rv, *t; |
3480 | 3621 |
3481 t = rv = rv_alloc(n); | 3622 t = rv = rv_alloc(n); |
3482 for(*t = *s++; *t; *t = *s++) t++; | 3623 for(*t = *s++; *t; *t = *s++) t++; |
3483 if (rve) | 3624 if (rve) |
3484 *rve = t; | 3625 *rve = t; |
3485 return rv; | 3626 return rv; |
3486 } | 3627 } |
(...skipping 102 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3589 spec_case, try_quick; | 3730 spec_case, try_quick; |
3590 Long L; | 3731 Long L; |
3591 #ifndef Sudden_Underflow | 3732 #ifndef Sudden_Underflow |
3592 int denorm; | 3733 int denorm; |
3593 ULong x; | 3734 ULong x; |
3594 #endif | 3735 #endif |
3595 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S; | 3736 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S; |
3596 U d2, eps, u; | 3737 U d2, eps, u; |
3597 double ds; | 3738 double ds; |
3598 char *s, *s0; | 3739 char *s, *s0; |
| 3740 #ifndef No_leftright |
| 3741 #ifdef IEEE_Arith |
| 3742 U eps1; |
| 3743 #endif |
| 3744 #endif |
3599 #ifdef SET_INEXACT | 3745 #ifdef SET_INEXACT |
3600 int inexact, oldinexact; | 3746 int inexact, oldinexact; |
3601 #endif | 3747 #endif |
3602 #ifdef Honor_FLT_ROUNDS /*{*/ | 3748 #ifdef Honor_FLT_ROUNDS /*{*/ |
3603 int Rounding; | 3749 int Rounding; |
3604 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ | 3750 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ |
3605 Rounding = Flt_Rounds; | 3751 Rounding = Flt_Rounds; |
3606 #else /*}{*/ | 3752 #else /*}{*/ |
3607 Rounding = 1; | 3753 Rounding = 1; |
3608 switch(fegetround()) { | 3754 switch(fegetround()) { |
(...skipping 52 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3661 if (Rounding >= 2) { | 3807 if (Rounding >= 2) { |
3662 if (*sign) | 3808 if (*sign) |
3663 Rounding = Rounding == 2 ? 0 : 2; | 3809 Rounding = Rounding == 2 ? 0 : 2; |
3664 else | 3810 else |
3665 if (Rounding != 2) | 3811 if (Rounding != 2) |
3666 Rounding = 0; | 3812 Rounding = 0; |
3667 } | 3813 } |
3668 #endif | 3814 #endif |
3669 | 3815 |
3670 b = d2b(&u, &be, &bbits); | 3816 b = d2b(&u, &be, &bbits); |
| 3817 #ifdef Sudden_Underflow |
3671 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); | 3818 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); |
3672 #ifndef Sudden_Underflow | 3819 #else |
3673 » if (i) { | 3820 » if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { |
3674 #endif | 3821 #endif |
3675 dval(&d2) = dval(&u); | 3822 dval(&d2) = dval(&u); |
3676 word0(&d2) &= Frac_mask1; | 3823 word0(&d2) &= Frac_mask1; |
3677 word0(&d2) |= Exp_11; | 3824 word0(&d2) |= Exp_11; |
3678 #ifdef IBM | 3825 #ifdef IBM |
3679 if (j = 11 - hi0bits(word0(&d2) & Frac_mask)) | 3826 if (j = 11 - hi0bits(word0(&d2) & Frac_mask)) |
3680 dval(&d2) /= 1 << j; | 3827 dval(&d2) /= 1 << j; |
3681 #endif | 3828 #endif |
3682 | 3829 |
3683 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 | 3830 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 |
(...skipping 134 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3818 dval(&u) /= bigtens[n_bigtens-1]; | 3965 dval(&u) /= bigtens[n_bigtens-1]; |
3819 ieps++; | 3966 ieps++; |
3820 } | 3967 } |
3821 for(; j; j >>= 1, i++) | 3968 for(; j; j >>= 1, i++) |
3822 if (j & 1) { | 3969 if (j & 1) { |
3823 ieps++; | 3970 ieps++; |
3824 ds *= bigtens[i]; | 3971 ds *= bigtens[i]; |
3825 } | 3972 } |
3826 dval(&u) /= ds; | 3973 dval(&u) /= ds; |
3827 } | 3974 } |
3828 » » else { | 3975 » » else if ((j1 = -k)) { |
3829 » » » j1 = -k; | 3976 » » » dval(&u) *= tens[j1 & 0xf]; |
3830 » » » if (j1) { | 3977 » » » for(j = j1 >> 4; j; j >>= 1, i++) |
3831 » » » » dval(&u) *= tens[j1 & 0xf]; | 3978 » » » » if (j & 1) { |
3832 » » » » for(j = j1 >> 4; j; j >>= 1, i++) | 3979 » » » » » ieps++; |
3833 » » » » » if (j & 1) { | 3980 » » » » » dval(&u) *= bigtens[i]; |
3834 » » » » » » ieps++; | 3981 » » » » » } |
3835 » » » » » » dval(&u) *= bigtens[i]; | |
3836 » » » » » » } | |
3837 » » » » } | |
3838 } | 3982 } |
3839 if (k_check && dval(&u) < 1. && ilim > 0) { | 3983 if (k_check && dval(&u) < 1. && ilim > 0) { |
3840 if (ilim1 <= 0) | 3984 if (ilim1 <= 0) |
3841 goto fast_failed; | 3985 goto fast_failed; |
3842 ilim = ilim1; | 3986 ilim = ilim1; |
3843 k--; | 3987 k--; |
3844 dval(&u) *= 10.; | 3988 dval(&u) *= 10.; |
3845 ieps++; | 3989 ieps++; |
3846 } | 3990 } |
3847 dval(&eps) = ieps*dval(&u) + 7.; | 3991 dval(&eps) = ieps*dval(&u) + 7.; |
3848 word0(&eps) -= (P-1)*Exp_msk1; | 3992 word0(&eps) -= (P-1)*Exp_msk1; |
3849 if (ilim == 0) { | 3993 if (ilim == 0) { |
3850 S = mhi = 0; | 3994 S = mhi = 0; |
3851 dval(&u) -= 5.; | 3995 dval(&u) -= 5.; |
3852 if (dval(&u) > dval(&eps)) | 3996 if (dval(&u) > dval(&eps)) |
3853 goto one_digit; | 3997 goto one_digit; |
3854 if (dval(&u) < -dval(&eps)) | 3998 if (dval(&u) < -dval(&eps)) |
3855 goto no_digits; | 3999 goto no_digits; |
3856 goto fast_failed; | 4000 goto fast_failed; |
3857 } | 4001 } |
3858 #ifndef No_leftright | 4002 #ifndef No_leftright |
3859 if (leftright) { | 4003 if (leftright) { |
3860 /* Use Steele & White method of only | 4004 /* Use Steele & White method of only |
3861 * generating digits needed. | 4005 * generating digits needed. |
3862 */ | 4006 */ |
3863 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); | 4007 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); |
| 4008 #ifdef IEEE_Arith |
| 4009 if (k0 < 0 && j1 >= 307) { |
| 4010 eps1.d = 1.01e256; /* 1.01 allows roundoff in th
e next few lines */ |
| 4011 word0(&eps1) -= Exp_msk1 * (Bias+P-1); |
| 4012 dval(&eps1) *= tens[j1 & 0xf]; |
| 4013 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++) |
| 4014 if (j & 1) |
| 4015 dval(&eps1) *= bigtens[i]; |
| 4016 if (eps.d < eps1.d) |
| 4017 eps.d = eps1.d; |
| 4018 } |
| 4019 #endif |
3864 for(i = 0;;) { | 4020 for(i = 0;;) { |
3865 » » » » L = (long)dval(&u); | 4021 » » » » L = dval(&u); |
3866 dval(&u) -= L; | 4022 dval(&u) -= L; |
3867 » » » » *s++ = '0' + (char)L; | 4023 » » » » *s++ = '0' + (int)L; |
| 4024 » » » » if (1. - dval(&u) < dval(&eps)) |
| 4025 » » » » » goto bump_up; |
3868 if (dval(&u) < dval(&eps)) | 4026 if (dval(&u) < dval(&eps)) |
3869 goto ret1; | 4027 goto ret1; |
3870 if (1. - dval(&u) < dval(&eps)) | |
3871 goto bump_up; | |
3872 if (++i >= ilim) | 4028 if (++i >= ilim) |
3873 break; | 4029 break; |
3874 dval(&eps) *= 10.; | 4030 dval(&eps) *= 10.; |
3875 dval(&u) *= 10.; | 4031 dval(&u) *= 10.; |
3876 } | 4032 } |
3877 } | 4033 } |
3878 else { | 4034 else { |
3879 #endif | 4035 #endif |
3880 /* Generate ilim digits, then fix them up. */ | 4036 /* Generate ilim digits, then fix them up. */ |
3881 dval(&eps) *= tens[ilim-1]; | 4037 dval(&eps) *= tens[ilim-1]; |
(...skipping 27 matching lines...) Expand all Loading... |
3909 | 4065 |
3910 if (be >= 0 && k <= Int_max) { | 4066 if (be >= 0 && k <= Int_max) { |
3911 /* Yes. */ | 4067 /* Yes. */ |
3912 ds = tens[k]; | 4068 ds = tens[k]; |
3913 if (ndigits < 0 && ilim <= 0) { | 4069 if (ndigits < 0 && ilim <= 0) { |
3914 S = mhi = 0; | 4070 S = mhi = 0; |
3915 if (ilim < 0 || dval(&u) <= 5*ds) | 4071 if (ilim < 0 || dval(&u) <= 5*ds) |
3916 goto no_digits; | 4072 goto no_digits; |
3917 goto one_digit; | 4073 goto one_digit; |
3918 } | 4074 } |
3919 » » for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) { | 4075 » » for(i = 1;; i++, dval(&u) *= 10.) { |
3920 L = (Long)(dval(&u) / ds); | 4076 L = (Long)(dval(&u) / ds); |
3921 dval(&u) -= L*ds; | 4077 dval(&u) -= L*ds; |
3922 #ifdef Check_FLT_ROUNDS | 4078 #ifdef Check_FLT_ROUNDS |
3923 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ | 4079 /* If FLT_ROUNDS == 2, L will usually be high by 1 */ |
3924 if (dval(&u) < 0) { | 4080 if (dval(&u) < 0) { |
3925 L--; | 4081 L--; |
3926 dval(&u) += ds; | 4082 dval(&u) += ds; |
3927 } | 4083 } |
3928 #endif | 4084 #endif |
3929 *s++ = '0' + (char)L; | 4085 *s++ = '0' + (char)L; |
3930 if (!dval(&u)) { | 4086 if (!dval(&u)) { |
3931 #ifdef SET_INEXACT | 4087 #ifdef SET_INEXACT |
3932 inexact = 0; | 4088 inexact = 0; |
3933 #endif | 4089 #endif |
3934 break; | 4090 break; |
3935 } | 4091 } |
3936 if (i == ilim) { | 4092 if (i == ilim) { |
3937 #ifdef Honor_FLT_ROUNDS | 4093 #ifdef Honor_FLT_ROUNDS |
3938 if (mode > 1) | 4094 if (mode > 1) |
3939 switch(Rounding) { | 4095 switch(Rounding) { |
3940 case 0: goto ret1; | 4096 case 0: goto ret1; |
3941 case 2: goto bump_up; | 4097 case 2: goto bump_up; |
3942 } | 4098 } |
3943 #endif | 4099 #endif |
3944 dval(&u) += dval(&u); | 4100 dval(&u) += dval(&u); |
3945 » » » » if (dval(&u) > ds || (dval(&u) == ds && L & 1))
{ | 4101 #ifdef ROUND_BIASED |
| 4102 » » » » if (dval(&u) >= ds) |
| 4103 #else |
| 4104 » » » » if (dval(&u) > ds || (dval(&u) == ds && L & 1)) |
| 4105 #endif |
| 4106 » » » » » { |
3946 bump_up: | 4107 bump_up: |
3947 while(*--s == '9') | 4108 while(*--s == '9') |
3948 if (s == s0) { | 4109 if (s == s0) { |
3949 k++; | 4110 k++; |
3950 *s = '0'; | 4111 *s = '0'; |
3951 break; | 4112 break; |
3952 } | 4113 } |
3953 ++*s++; | 4114 ++*s++; |
3954 } | 4115 } |
3955 break; | 4116 break; |
(...skipping 64 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
4020 } | 4181 } |
4021 } | 4182 } |
4022 | 4183 |
4023 /* Arrange for convenient computation of quotients: | 4184 /* Arrange for convenient computation of quotients: |
4024 * shift left if necessary so divisor has 4 leading 0 bits. | 4185 * shift left if necessary so divisor has 4 leading 0 bits. |
4025 * | 4186 * |
4026 * Perhaps we should just compute leading 28 bits of S once | 4187 * Perhaps we should just compute leading 28 bits of S once |
4027 * and for all and pass them and a shift to quorem, so it | 4188 * and for all and pass them and a shift to quorem, so it |
4028 * can do shifts and ors to compute the numerator for q. | 4189 * can do shifts and ors to compute the numerator for q. |
4029 */ | 4190 */ |
4030 #ifdef Pack_32 | |
4031 i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f; | |
4032 if (i) | |
4033 i = 32 - i; | |
4034 #define iInc 28 | |
4035 #else | |
4036 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) | |
4037 i = 16 - i; | |
4038 #define iInc 12 | |
4039 #endif | |
4040 i = dshift(S, s2); | 4191 i = dshift(S, s2); |
4041 b2 += i; | 4192 b2 += i; |
4042 m2 += i; | 4193 m2 += i; |
4043 s2 += i; | 4194 s2 += i; |
4044 if (b2 > 0) | 4195 if (b2 > 0) |
4045 b = lshift(b, b2); | 4196 b = lshift(b, b2); |
4046 if (s2 > 0) | 4197 if (s2 > 0) |
4047 S = lshift(S, s2); | 4198 S = lshift(S, s2); |
4048 if (k_check) { | 4199 if (k_check) { |
4049 if (cmp(b,S) < 0) { | 4200 if (cmp(b,S) < 0) { |
(...skipping 72 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
4122 #ifdef Honor_FLT_ROUNDS | 4273 #ifdef Honor_FLT_ROUNDS |
4123 if (mode > 1) | 4274 if (mode > 1) |
4124 switch(Rounding) { | 4275 switch(Rounding) { |
4125 case 0: goto accept_dig; | 4276 case 0: goto accept_dig; |
4126 case 2: goto keep_dig; | 4277 case 2: goto keep_dig; |
4127 } | 4278 } |
4128 #endif /*Honor_FLT_ROUNDS*/ | 4279 #endif /*Honor_FLT_ROUNDS*/ |
4129 if (j1 > 0) { | 4280 if (j1 > 0) { |
4130 b = lshift(b, 1); | 4281 b = lshift(b, 1); |
4131 j1 = cmp(b, S); | 4282 j1 = cmp(b, S); |
| 4283 #ifdef ROUND_BIASED |
| 4284 if (j1 >= 0 /*)*/ |
| 4285 #else |
4132 if ((j1 > 0 || (j1 == 0 && dig & 1)) | 4286 if ((j1 > 0 || (j1 == 0 && dig & 1)) |
| 4287 #endif |
4133 && dig++ == '9') | 4288 && dig++ == '9') |
4134 goto round_9_up; | 4289 goto round_9_up; |
4135 } | 4290 } |
4136 accept_dig: | 4291 accept_dig: |
4137 *s++ = (char)dig; | 4292 *s++ = (char)dig; |
4138 goto ret; | 4293 goto ret; |
4139 } | 4294 } |
4140 if (j1 > 0) { | 4295 if (j1 > 0) { |
4141 #ifdef Honor_FLT_ROUNDS | 4296 #ifdef Honor_FLT_ROUNDS |
4142 if (!Rounding) | 4297 if (!Rounding) |
(...skipping 40 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
4183 /* Round off last digit */ | 4338 /* Round off last digit */ |
4184 | 4339 |
4185 #ifdef Honor_FLT_ROUNDS | 4340 #ifdef Honor_FLT_ROUNDS |
4186 switch(Rounding) { | 4341 switch(Rounding) { |
4187 case 0: goto trimzeros; | 4342 case 0: goto trimzeros; |
4188 case 2: goto roundoff; | 4343 case 2: goto roundoff; |
4189 } | 4344 } |
4190 #endif | 4345 #endif |
4191 b = lshift(b, 1); | 4346 b = lshift(b, 1); |
4192 j = cmp(b, S); | 4347 j = cmp(b, S); |
4193 » if (j > 0 || (j == 0 && dig & 1)) { | 4348 #ifdef ROUND_BIASED |
| 4349 » if (j >= 0) |
| 4350 #else |
| 4351 » if (j > 0 || (j == 0 && dig & 1)) |
| 4352 #endif |
| 4353 » » { |
4194 roundoff: | 4354 roundoff: |
4195 while(*--s == '9') | 4355 while(*--s == '9') |
4196 if (s == s0) { | 4356 if (s == s0) { |
4197 k++; | 4357 k++; |
4198 *s++ = '1'; | 4358 *s++ = '1'; |
4199 goto ret; | 4359 goto ret; |
4200 } | 4360 } |
4201 ++*s++; | 4361 ++*s++; |
4202 } | 4362 } |
4203 else { | 4363 else { |
(...skipping 24 matching lines...) Expand all Loading... |
4228 #endif | 4388 #endif |
4229 Bfree(b); | 4389 Bfree(b); |
4230 *s = 0; | 4390 *s = 0; |
4231 *decpt = k + 1; | 4391 *decpt = k + 1; |
4232 if (rve) | 4392 if (rve) |
4233 *rve = s; | 4393 *rve = s; |
4234 return s0; | 4394 return s0; |
4235 } | 4395 } |
4236 | 4396 |
4237 } // namespace dmg_fp | 4397 } // namespace dmg_fp |
OLD | NEW |