1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998-2001 by Lucent Technologies
8 Permission to use, copy, modify, and distribute this software and
9 its documentation for any purpose and without fee is hereby
10 granted, provided that the above copyright notice appear in all
11 copies and that both that the copyright notice and this
12 permission notice and warranty disclaimer appear in supporting
13 documentation, and that the name of Lucent or any of its entities
14 not be used in advertising or publicity pertaining to
15 distribution of the software without specific, written prior
18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
43 #define Avoid_Underflow
45 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
46 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
47 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
48 9007199254740992.*9007199254740992.e-256
53 #ifdef Honor_FLT_ROUNDS
54 #undef Check_FLT_ROUNDS
55 #define Check_FLT_ROUNDS
57 #define Rounding Flt_Rounds
60 #ifdef Avoid_Underflow /*{*/
64 (x, scale) U *x; int scale;
74 if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
75 return rv; /* Is there an example where i <= 0 ? */
76 word0(&u) = Exp_1 + (i << Exp_shift);
85 (s00, se) CONST char *s00; char **se;
87 (CONST char *s00, char **se)
90 #ifdef Avoid_Underflow
93 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
94 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
95 CONST char *s, *s0, *s1;
98 U adj, aadj1, rv, rv0;
100 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
101 #ifdef Avoid_Underflow
105 int inexact, oldinexact;
107 #ifdef USE_LOCALE /*{{*/
108 #ifdef NO_LOCALE_CACHE
109 char *decimalpoint = localeconv()->decimal_point;
110 int dplen = strlen(decimalpoint);
113 static char *decimalpoint_cache;
115 if (!(s0 = decimalpoint_cache)) {
116 s0 = localeconv()->decimal_point;
117 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
118 strcpy(decimalpoint_cache, s0);
119 s0 = decimalpoint_cache;
123 decimalpoint = (char*)s0;
124 #endif /*NO_LOCALE_CACHE*/
125 #else /*USE_LOCALE}{*/
127 #endif /*USE_LOCALE}}*/
129 #ifdef Honor_FLT_ROUNDS /*{*/
131 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
132 Rounding = Flt_Rounds;
135 switch(fegetround()) {
136 case FE_TOWARDZERO: Rounding = 0; break;
137 case FE_UPWARD: Rounding = 2; break;
138 case FE_DOWNWARD: Rounding = 3;
143 sign = nz0 = nz = decpt = 0;
145 for(s = s00;;s++) switch(*s) {
167 #ifndef NO_HEX_FP /*{*/
169 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
176 #ifdef Honor_FLT_ROUNDS
178 fpi1.rounding = Rounding;
182 switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
183 case STRTOG_NoNumber:
190 copybits(bits, fpi.nbits, bb);
193 ULtod(((U*)&rv)->L, bits, exp, i);
206 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
213 if (c == *decimalpoint) {
214 for(i = 1; decimalpoint[i]; ++i)
215 if (s[i] != decimalpoint[i])
225 for(; c == '0'; c = *++s)
227 if (c > '0' && c <= '9') {
235 for(; c >= '0' && c <= '9'; c = *++s) {
240 for(i = 1; i < nz; i++)
243 else if (nd <= DBL_DIG + 1)
247 else if (nd <= DBL_DIG + 1)
255 if (c == 'e' || c == 'E') {
256 if (!nd && !nz && !nz0) {
267 if (c >= '0' && c <= '9') {
270 if (c > '0' && c <= '9') {
273 while((c = *++s) >= '0' && c <= '9')
275 if (s - s1 > 8 || L > 19999)
276 /* Avoid confusion from exponents
277 * so large that e might overflow.
279 e = 19999; /* safe for 16 bit ints */
294 /* Check for Nan and Infinity */
296 static FPI fpinan = /* only 52 explicit bits */
297 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
302 if (match(&s,"nf")) {
304 if (!match(&s,"inity"))
306 word0(&rv) = 0x7ff00000;
313 if (match(&s, "an")) {
316 && hexnan(&s, &fpinan, bits)
318 word0(&rv) = 0x7ff00000 | bits[1];
319 word1(&rv) = bits[0];
323 word0(&rv) = NAN_WORD0;
324 word1(&rv) = NAN_WORD1;
331 #endif /* INFNAN_CHECK */
340 /* Now we have nd0 digits, starting at s0, followed by a
341 * decimal point, followed by nd-nd0 digits. The number we're
342 * after is the integer represented by those digits times
347 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
352 oldinexact = get_inexact();
354 dval(&rv) = tens[k - 9] * dval(&rv) + z;
359 #ifndef Honor_FLT_ROUNDS
371 #ifdef Honor_FLT_ROUNDS
372 /* round correctly FLT_ROUNDS = 2 or 3 */
378 /* rv = */ rounded_product(dval(&rv), tens[e]);
383 if (e <= Ten_pmax + i) {
384 /* A fancier test would sometimes let us do
385 * this for larger i values.
387 #ifdef Honor_FLT_ROUNDS
388 /* round correctly FLT_ROUNDS = 2 or 3 */
395 dval(&rv) *= tens[i];
397 /* VAX exponent range is so narrow we must
398 * worry about overflow here...
401 word0(&rv) -= P*Exp_msk1;
402 /* rv = */ rounded_product(dval(&rv), tens[e]);
403 if ((word0(&rv) & Exp_mask)
404 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
406 word0(&rv) += P*Exp_msk1;
408 /* rv = */ rounded_product(dval(&rv), tens[e]);
413 #ifndef Inaccurate_Divide
414 else if (e >= -Ten_pmax) {
415 #ifdef Honor_FLT_ROUNDS
416 /* round correctly FLT_ROUNDS = 2 or 3 */
422 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
433 oldinexact = get_inexact();
435 #ifdef Avoid_Underflow
438 #ifdef Honor_FLT_ROUNDS
441 Rounding = Rounding == 2 ? 0 : 2;
447 #endif /*IEEE_Arith*/
449 /* Get starting approximation = rv * 10**e1 */
452 if ( (i = e1 & 15) !=0)
453 dval(&rv) *= tens[i];
455 if (e1 > DBL_MAX_10_EXP) {
457 /* Can't trust HUGE_VAL */
459 #ifdef Honor_FLT_ROUNDS
461 case 0: /* toward 0 */
462 case 3: /* toward -infinity */
467 word0(&rv) = Exp_mask;
470 #else /*Honor_FLT_ROUNDS*/
471 word0(&rv) = Exp_mask;
473 #endif /*Honor_FLT_ROUNDS*/
475 /* set overflow bit */
477 dval(&rv0) *= dval(&rv0);
482 #endif /*IEEE_Arith*/
497 for(j = 0; e1 > 1; j++, e1 >>= 1)
499 dval(&rv) *= bigtens[j];
500 /* The last multiplication could overflow. */
501 word0(&rv) -= P*Exp_msk1;
502 dval(&rv) *= bigtens[j];
503 if ((z = word0(&rv) & Exp_mask)
504 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
506 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
507 /* set to largest number */
508 /* (Can't trust DBL_MAX) */
513 word0(&rv) += P*Exp_msk1;
518 if ( (i = e1 & 15) !=0)
519 dval(&rv) /= tens[i];
521 if (e1 >= 1 << n_bigtens)
523 #ifdef Avoid_Underflow
526 for(j = 0; e1 > 0; j++, e1 >>= 1)
528 dval(&rv) *= tinytens[j];
529 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
530 >> Exp_shift)) > 0) {
531 /* scaled rv is denormal; zap j low bits */
535 word0(&rv) = (P+2)*Exp_msk1;
537 word0(&rv) &= 0xffffffff << (j-32);
540 word1(&rv) &= 0xffffffff << j;
543 for(j = 0; e1 > 1; j++, e1 >>= 1)
545 dval(&rv) *= tinytens[j];
546 /* The last multiplication could underflow. */
547 dval(&rv0) = dval(&rv);
548 dval(&rv) *= tinytens[j];
550 dval(&rv) = 2.*dval(&rv0);
551 dval(&rv) *= tinytens[j];
558 #ifndef Avoid_Underflow
561 /* The refinement below will clean
562 * this approximation up.
569 /* Now the hard part -- adjusting rv to the correct value.*/
571 /* Put digits into bd: true value = bd * 10^e */
573 bd0 = s2b(s0, nd0, nd, y, dplen);
578 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
594 #ifdef Honor_FLT_ROUNDS
598 #ifdef Avoid_Underflow
602 i = j + bbbits - 1; /* logb(rv) */
604 if (i < Emin) { /* denormal */
610 Lsb1 = Lsb << (i-32);
612 #else /*Avoid_Underflow*/
613 #ifdef Sudden_Underflow
615 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
619 #else /*Sudden_Underflow*/
621 i = j + bbbits - 1; /* logb(&rv) */
622 if (i < Emin) /* denormal */
626 #endif /*Sudden_Underflow*/
627 #endif /*Avoid_Underflow*/
630 #ifdef Avoid_Underflow
633 i = bb2 < bd2 ? bb2 : bd2;
642 bs = pow5mult(bs, bb5);
648 bb = lshift(bb, bb2);
650 bd = pow5mult(bd, bd5);
652 bd = lshift(bd, bd2);
654 bs = lshift(bs, bs2);
655 delta = diff(bb, bd);
659 #ifdef Honor_FLT_ROUNDS
662 /* Error is less than an ulp */
663 if (!delta->x[0] && delta->wds <= 1) {
679 && !(word0(&rv) & Frac_mask)) {
680 y = word0(&rv) & Exp_mask;
681 #ifdef Avoid_Underflow
682 if (!scale || y > 2*P*Exp_msk1)
687 delta = lshift(delta,Log2P);
688 if (cmp(delta, bs) <= 0)
693 #ifdef Avoid_Underflow
694 if (scale && (y = word0(&rv) & Exp_mask)
696 word0(&adj) += (2*P+1)*Exp_msk1 - y;
698 #ifdef Sudden_Underflow
699 if ((word0(&rv) & Exp_mask) <=
701 word0(&rv) += P*Exp_msk1;
702 dval(&rv) += adj*ulp(&rv);
703 word0(&rv) -= P*Exp_msk1;
706 #endif /*Sudden_Underflow*/
707 #endif /*Avoid_Underflow*/
708 dval(&rv) += adj.d*ulp(&rv);
712 dval(&adj) = ratio(delta, bs);
715 if (adj.d <= 0x7ffffffe) {
716 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
719 if (!((Rounding>>1) ^ dsign))
724 #ifdef Avoid_Underflow
725 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
726 word0(&adj) += (2*P+1)*Exp_msk1 - y;
728 #ifdef Sudden_Underflow
729 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
730 word0(&rv) += P*Exp_msk1;
731 dval(&adj) *= ulp(&rv);
736 word0(&rv) -= P*Exp_msk1;
739 #endif /*Sudden_Underflow*/
740 #endif /*Avoid_Underflow*/
741 dval(&adj) *= ulp(&rv);
743 if (word0(&rv) == Big0 && word1(&rv) == Big1)
751 #endif /*Honor_FLT_ROUNDS*/
754 /* Error is less than half an ulp -- check for
755 * special case of mantissa a power of two.
757 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
759 #ifdef Avoid_Underflow
760 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
762 || (word0(&rv) & Exp_mask) <= Exp_msk1
767 if (!delta->x[0] && delta->wds <= 1)
772 if (!delta->x[0] && delta->wds <= 1) {
779 delta = lshift(delta,Log2P);
780 if (cmp(delta, bs) > 0)
785 /* exactly half-way between */
787 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
789 #ifdef Avoid_Underflow
790 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
791 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
794 /*boundary case -- increment exponent*/
795 if (word0(&rv) == Big0 && word1(&rv) == Big1)
797 word0(&rv) = (word0(&rv) & Exp_mask)
804 #ifdef Avoid_Underflow
810 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
812 /* boundary case -- decrement exponent */
813 #ifdef Sudden_Underflow /*{{*/
814 L = word0(&rv) & Exp_mask;
818 #ifdef Avoid_Underflow
819 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
822 #endif /*Avoid_Underflow*/
826 #else /*Sudden_Underflow}{*/
827 #ifdef Avoid_Underflow
829 L = word0(&rv) & Exp_mask;
830 if (L <= (2*P+1)*Exp_msk1) {
831 if (L > (P+2)*Exp_msk1)
835 /* rv = smallest denormal */
839 #endif /*Avoid_Underflow*/
840 L = (word0(&rv) & Exp_mask) - Exp_msk1;
841 #endif /*Sudden_Underflow}}*/
842 word0(&rv) = L | Bndry_mask1;
843 word1(&rv) = 0xffffffff;
851 #ifdef Avoid_Underflow
853 if (!(word0(&rv) & Lsb1))
856 else if (!(word1(&rv) & Lsb))
859 if (!(word1(&rv) & LSB))
864 #ifdef Avoid_Underflow
865 dval(&rv) += sulp(&rv, scale);
867 dval(&rv) += ulp(&rv);
871 #ifdef Avoid_Underflow
872 dval(&rv) -= sulp(&rv, scale);
874 dval(&rv) -= ulp(&rv);
876 #ifndef Sudden_Underflow
881 #ifdef Avoid_Underflow
887 if ((aadj = ratio(delta, bs)) <= 2.) {
889 aadj = dval(&aadj1) = 1.;
890 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
891 #ifndef Sudden_Underflow
892 if (word1(&rv) == Tiny1 && !word0(&rv))
899 /* special case -- power of FLT_RADIX to be */
900 /* rounded down... */
902 if (aadj < 2./FLT_RADIX)
906 dval(&aadj1) = -aadj;
911 dval(&aadj1) = dsign ? aadj : -aadj;
912 #ifdef Check_FLT_ROUNDS
914 case 2: /* towards +infinity */
917 case 0: /* towards 0 */
918 case 3: /* towards -infinity */
924 #endif /*Check_FLT_ROUNDS*/
926 y = word0(&rv) & Exp_mask;
928 /* Check for overflow */
930 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
931 dval(&rv0) = dval(&rv);
932 word0(&rv) -= P*Exp_msk1;
933 dval(&adj) = dval(&aadj1) * ulp(&rv);
934 dval(&rv) += dval(&adj);
935 if ((word0(&rv) & Exp_mask) >=
936 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
937 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
944 word0(&rv) += P*Exp_msk1;
947 #ifdef Avoid_Underflow
948 if (scale && y <= 2*P*Exp_msk1) {
949 if (aadj <= 0x7fffffff) {
953 dval(&aadj1) = dsign ? aadj : -aadj;
955 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
957 dval(&adj) = dval(&aadj1) * ulp(&rv);
958 dval(&rv) += dval(&adj);
960 #ifdef Sudden_Underflow
961 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
962 dval(&rv0) = dval(&rv);
963 word0(&rv) += P*Exp_msk1;
964 dval(&adj) = dval(&aadj1) * ulp(&rv);
967 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
969 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
972 if (word0(&rv0) == Tiny0
973 && word1(&rv0) == Tiny1)
980 word0(&rv) -= P*Exp_msk1;
983 dval(&adj) = dval(&aadj1) * ulp(&rv);
986 #else /*Sudden_Underflow*/
987 /* Compute dval(&adj) so that the IEEE rounding rules will
988 * correctly round rv + dval(&adj) in some half-way cases.
989 * If rv * ulp(&rv) is denormalized (i.e.,
990 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
991 * trouble from bits lost to denormalization;
992 * example: 1.2e-307 .
994 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
995 dval(&aadj1) = (double)(int)(aadj + 0.5);
997 dval(&aadj1) = -dval(&aadj1);
999 dval(&adj) = dval(&aadj1) * ulp(&rv);
1001 #endif /*Sudden_Underflow*/
1002 #endif /*Avoid_Underflow*/
1004 z = word0(&rv) & Exp_mask;
1006 #ifdef Avoid_Underflow
1010 /* Can we stop now? */
1013 /* The tolerances below are conservative. */
1014 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1015 if (aadj < .4999999 || aadj > .5000001)
1018 else if (aadj < .4999999/FLT_RADIX)
1036 word0(&rv0) = Exp_1 + (70 << Exp_shift);
1041 else if (!oldinexact)
1044 #ifdef Avoid_Underflow
1046 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1048 dval(&rv) *= dval(&rv0);
1050 /* try to avoid the bug of testing an 8087 register value */
1052 if (!(word0(&rv) & Exp_mask))
1054 if (word0(&rv) == 0 && word1(&rv) == 0)
1059 #endif /* Avoid_Underflow */
1061 if (inexact && !(word0(&rv) & Exp_mask)) {
1062 /* set underflow bit */
1063 dval(&rv0) = 1e-300;
1064 dval(&rv0) *= dval(&rv0);
1070 return sign ? -dval(&rv) : dval(&rv);