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 "."). */
45 #define Avoid_Underflow
47 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
48 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
49 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
50 9007199254740992.*9007199254740992.e-256
55 #ifdef Honor_FLT_ROUNDS
56 #undef Check_FLT_ROUNDS
57 #define Check_FLT_ROUNDS
59 #define Rounding Flt_Rounds
65 (s00, se) CONST char *s00; char **se;
67 (CONST char *s00, char **se)
70 #ifdef Avoid_Underflow
73 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
74 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
75 CONST char *s, *s0, *s1;
76 double aadj, aadj1, adj, rv, rv0;
79 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
81 int inexact, oldinexact;
83 #ifdef USE_LOCALE /*{{*/
84 #ifdef NO_LOCALE_CACHE
85 char *decimalpoint = localeconv()->decimal_point;
86 int dplen = strlen(decimalpoint);
89 static char *decimalpoint_cache;
91 if (!(s0 = decimalpoint_cache)) {
92 s0 = localeconv()->decimal_point;
93 if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
94 strcpy(decimalpoint_cache, s0);
95 s0 = decimalpoint_cache;
99 decimalpoint = (char*)s0;
100 #endif /*NO_LOCALE_CACHE*/
101 #else /*USE_LOCALE}{*/
103 #endif /*USE_LOCALE}}*/
105 #ifdef Honor_FLT_ROUNDS /*{*/
107 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
108 Rounding = Flt_Rounds;
111 switch(fegetround()) {
112 case FE_TOWARDZERO: Rounding = 0; break;
113 case FE_UPWARD: Rounding = 2; break;
114 case FE_DOWNWARD: Rounding = 3;
119 sign = nz0 = nz = decpt = 0;
121 for(s = s00;;s++) switch(*s) {
143 #ifndef NO_HEX_FP /*{*/
145 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
152 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/
154 #ifdef Honor_FLT_ROUNDS /*{{*/
155 fpi1.rounding = Rounding;
157 switch(fegetround()) {
158 case FE_TOWARDZERO: fpi1.rounding = 0; break;
159 case FE_UPWARD: fpi1.rounding = 2; break;
160 case FE_DOWNWARD: fpi1.rounding = 3;
166 switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
167 case STRTOG_NoNumber:
174 copybits(bits, fpi.nbits, bb);
177 ULtod(((U*)&rv)->L, bits, exp, i);
190 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
197 if (c == *decimalpoint) {
198 for(i = 1; decimalpoint[i]; ++i)
199 if (s[i] != decimalpoint[i])
209 for(; c == '0'; c = *++s)
211 if (c > '0' && c <= '9') {
219 for(; c >= '0' && c <= '9'; c = *++s) {
224 for(i = 1; i < nz; i++)
227 else if (nd <= DBL_DIG + 1)
231 else if (nd <= DBL_DIG + 1)
239 if (c == 'e' || c == 'E') {
240 if (!nd && !nz && !nz0) {
251 if (c >= '0' && c <= '9') {
254 if (c > '0' && c <= '9') {
257 while((c = *++s) >= '0' && c <= '9')
259 if (s - s1 > 8 || L > 19999)
260 /* Avoid confusion from exponents
261 * so large that e might overflow.
263 e = 19999; /* safe for 16 bit ints */
278 /* Check for Nan and Infinity */
280 static FPI fpinan = /* only 52 explicit bits */
281 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
286 if (match(&s,"nf")) {
288 if (!match(&s,"inity"))
290 word0(rv) = 0x7ff00000;
297 if (match(&s, "an")) {
300 && hexnan(&s, &fpinan, bits)
302 word0(rv) = 0x7ff80000 | bits[1];
307 word0(rv) = NAN_WORD0;
308 word1(rv) = NAN_WORD1;
315 #endif /* INFNAN_CHECK */
324 /* Now we have nd0 digits, starting at s0, followed by a
325 * decimal point, followed by nd-nd0 digits. The number we're
326 * after is the integer represented by those digits times
331 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
336 oldinexact = get_inexact();
338 dval(rv) = tens[k - 9] * dval(rv) + z;
343 #ifndef Honor_FLT_ROUNDS
355 #ifdef Honor_FLT_ROUNDS
356 /* round correctly FLT_ROUNDS = 2 or 3 */
362 /* rv = */ rounded_product(dval(rv), tens[e]);
367 if (e <= Ten_pmax + i) {
368 /* A fancier test would sometimes let us do
369 * this for larger i values.
371 #ifdef Honor_FLT_ROUNDS
372 /* round correctly FLT_ROUNDS = 2 or 3 */
381 /* VAX exponent range is so narrow we must
382 * worry about overflow here...
385 word0(rv) -= P*Exp_msk1;
386 /* rv = */ rounded_product(dval(rv), tens[e]);
387 if ((word0(rv) & Exp_mask)
388 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
390 word0(rv) += P*Exp_msk1;
392 /* rv = */ rounded_product(dval(rv), tens[e]);
397 #ifndef Inaccurate_Divide
398 else if (e >= -Ten_pmax) {
399 #ifdef Honor_FLT_ROUNDS
400 /* round correctly FLT_ROUNDS = 2 or 3 */
406 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
417 oldinexact = get_inexact();
419 #ifdef Avoid_Underflow
422 #ifdef Honor_FLT_ROUNDS
425 Rounding = Rounding == 2 ? 0 : 2;
431 #endif /*IEEE_Arith*/
433 /* Get starting approximation = rv * 10**e1 */
436 if ( (i = e1 & 15) !=0)
439 if (e1 > DBL_MAX_10_EXP) {
444 /* Can't trust HUGE_VAL */
446 #ifdef Honor_FLT_ROUNDS
448 case 0: /* toward 0 */
449 case 3: /* toward -infinity */
454 word0(rv) = Exp_mask;
457 #else /*Honor_FLT_ROUNDS*/
458 word0(rv) = Exp_mask;
460 #endif /*Honor_FLT_ROUNDS*/
462 /* set overflow bit */
464 dval(rv0) *= dval(rv0);
469 #endif /*IEEE_Arith*/
475 for(j = 0; e1 > 1; j++, e1 >>= 1)
477 dval(rv) *= bigtens[j];
478 /* The last multiplication could overflow. */
479 word0(rv) -= P*Exp_msk1;
480 dval(rv) *= bigtens[j];
481 if ((z = word0(rv) & Exp_mask)
482 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
484 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
485 /* set to largest number */
486 /* (Can't trust DBL_MAX) */
491 word0(rv) += P*Exp_msk1;
496 if ( (i = e1 & 15) !=0)
499 if (e1 >= 1 << n_bigtens)
501 #ifdef Avoid_Underflow
504 for(j = 0; e1 > 0; j++, e1 >>= 1)
506 dval(rv) *= tinytens[j];
507 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
508 >> Exp_shift)) > 0) {
509 /* scaled rv is denormal; zap j low bits */
513 word0(rv) = (P+2)*Exp_msk1;
515 word0(rv) &= 0xffffffff << j-32;
518 word1(rv) &= 0xffffffff << j;
521 for(j = 0; e1 > 1; j++, e1 >>= 1)
523 dval(rv) *= tinytens[j];
524 /* The last multiplication could underflow. */
525 dval(rv0) = dval(rv);
526 dval(rv) *= tinytens[j];
528 dval(rv) = 2.*dval(rv0);
529 dval(rv) *= tinytens[j];
541 #ifndef Avoid_Underflow
544 /* The refinement below will clean
545 * this approximation up.
552 /* Now the hard part -- adjusting rv to the correct value.*/
554 /* Put digits into bd: true value = bd * 10^e */
556 bd0 = s2b(s0, nd0, nd, y, dplen);
561 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
577 #ifdef Honor_FLT_ROUNDS
581 #ifdef Avoid_Underflow
583 i = j + bbbits - 1; /* logb(rv) */
584 if (i < Emin) /* denormal */
588 #else /*Avoid_Underflow*/
589 #ifdef Sudden_Underflow
591 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
595 #else /*Sudden_Underflow*/
597 i = j + bbbits - 1; /* logb(rv) */
598 if (i < Emin) /* denormal */
602 #endif /*Sudden_Underflow*/
603 #endif /*Avoid_Underflow*/
606 #ifdef Avoid_Underflow
609 i = bb2 < bd2 ? bb2 : bd2;
618 bs = pow5mult(bs, bb5);
624 bb = lshift(bb, bb2);
626 bd = pow5mult(bd, bd5);
628 bd = lshift(bd, bd2);
630 bs = lshift(bs, bs2);
631 delta = diff(bb, bd);
635 #ifdef Honor_FLT_ROUNDS
638 /* Error is less than an ulp */
639 if (!delta->x[0] && delta->wds <= 1) {
655 && !(word0(rv) & Frac_mask)) {
656 y = word0(rv) & Exp_mask;
657 #ifdef Avoid_Underflow
658 if (!scale || y > 2*P*Exp_msk1)
663 delta = lshift(delta,Log2P);
664 if (cmp(delta, bs) <= 0)
669 #ifdef Avoid_Underflow
670 if (scale && (y = word0(rv) & Exp_mask)
672 word0(adj) += (2*P+1)*Exp_msk1 - y;
674 #ifdef Sudden_Underflow
675 if ((word0(rv) & Exp_mask) <=
677 word0(rv) += P*Exp_msk1;
678 dval(rv) += adj*ulp(dval(rv));
679 word0(rv) -= P*Exp_msk1;
682 #endif /*Sudden_Underflow*/
683 #endif /*Avoid_Underflow*/
684 dval(rv) += adj*ulp(dval(rv));
688 adj = ratio(delta, bs);
691 if (adj <= 0x7ffffffe) {
692 /* adj = Rounding ? ceil(adj) : floor(adj); */
695 if (!((Rounding>>1) ^ dsign))
700 #ifdef Avoid_Underflow
701 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
702 word0(adj) += (2*P+1)*Exp_msk1 - y;
704 #ifdef Sudden_Underflow
705 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
706 word0(rv) += P*Exp_msk1;
707 adj *= ulp(dval(rv));
712 word0(rv) -= P*Exp_msk1;
715 #endif /*Sudden_Underflow*/
716 #endif /*Avoid_Underflow*/
717 adj *= ulp(dval(rv));
719 if (word0(rv) == Big0 && word1(rv) == Big1)
727 #endif /*Honor_FLT_ROUNDS*/
730 /* Error is less than half an ulp -- check for
731 * special case of mantissa a power of two.
733 if (dsign || word1(rv) || word0(rv) & Bndry_mask
735 #ifdef Avoid_Underflow
736 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
738 || (word0(rv) & Exp_mask) <= Exp_msk1
743 if (!delta->x[0] && delta->wds <= 1)
748 if (!delta->x[0] && delta->wds <= 1) {
755 delta = lshift(delta,Log2P);
756 if (cmp(delta, bs) > 0)
761 /* exactly half-way between */
763 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
765 #ifdef Avoid_Underflow
766 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
767 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
770 /*boundary case -- increment exponent*/
771 word0(rv) = (word0(rv) & Exp_mask)
778 #ifdef Avoid_Underflow
784 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
786 /* boundary case -- decrement exponent */
787 #ifdef Sudden_Underflow /*{{*/
788 L = word0(rv) & Exp_mask;
792 #ifdef Avoid_Underflow
793 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
796 #endif /*Avoid_Underflow*/
800 #else /*Sudden_Underflow}{*/
801 #ifdef Avoid_Underflow
803 L = word0(rv) & Exp_mask;
804 if (L <= (2*P+1)*Exp_msk1) {
805 if (L > (P+2)*Exp_msk1)
809 /* rv = smallest denormal */
813 #endif /*Avoid_Underflow*/
814 L = (word0(rv) & Exp_mask) - Exp_msk1;
815 #endif /*Sudden_Underflow}}*/
816 word0(rv) = L | Bndry_mask1;
817 word1(rv) = 0xffffffff;
825 if (!(word1(rv) & LSB))
829 dval(rv) += ulp(dval(rv));
832 dval(rv) -= ulp(dval(rv));
833 #ifndef Sudden_Underflow
838 #ifdef Avoid_Underflow
844 if ((aadj = ratio(delta, bs)) <= 2.) {
847 else if (word1(rv) || word0(rv) & Bndry_mask) {
848 #ifndef Sudden_Underflow
849 if (word1(rv) == Tiny1 && !word0(rv))
856 /* special case -- power of FLT_RADIX to be */
857 /* rounded down... */
859 if (aadj < 2./FLT_RADIX)
868 aadj1 = dsign ? aadj : -aadj;
869 #ifdef Check_FLT_ROUNDS
871 case 2: /* towards +infinity */
874 case 0: /* towards 0 */
875 case 3: /* towards -infinity */
881 #endif /*Check_FLT_ROUNDS*/
883 y = word0(rv) & Exp_mask;
885 /* Check for overflow */
887 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
888 dval(rv0) = dval(rv);
889 word0(rv) -= P*Exp_msk1;
890 adj = aadj1 * ulp(dval(rv));
892 if ((word0(rv) & Exp_mask) >=
893 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
894 if (word0(rv0) == Big0 && word1(rv0) == Big1)
901 word0(rv) += P*Exp_msk1;
904 #ifdef Avoid_Underflow
905 if (scale && y <= 2*P*Exp_msk1) {
906 if (aadj <= 0x7fffffff) {
910 aadj1 = dsign ? aadj : -aadj;
912 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
914 adj = aadj1 * ulp(dval(rv));
917 #ifdef Sudden_Underflow
918 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
919 dval(rv0) = dval(rv);
920 word0(rv) += P*Exp_msk1;
921 adj = aadj1 * ulp(dval(rv));
924 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
926 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
929 if (word0(rv0) == Tiny0
930 && word1(rv0) == Tiny1)
937 word0(rv) -= P*Exp_msk1;
940 adj = aadj1 * ulp(dval(rv));
943 #else /*Sudden_Underflow*/
944 /* Compute adj so that the IEEE rounding rules will
945 * correctly round rv + adj in some half-way cases.
946 * If rv * ulp(rv) is denormalized (i.e.,
947 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
948 * trouble from bits lost to denormalization;
949 * example: 1.2e-307 .
951 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
952 aadj1 = (double)(int)(aadj + 0.5);
956 adj = aadj1 * ulp(dval(rv));
958 #endif /*Sudden_Underflow*/
959 #endif /*Avoid_Underflow*/
961 z = word0(rv) & Exp_mask;
963 #ifdef Avoid_Underflow
967 /* Can we stop now? */
970 /* The tolerances below are conservative. */
971 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
972 if (aadj < .4999999 || aadj > .5000001)
975 else if (aadj < .4999999/FLT_RADIX)
988 word0(rv0) = Exp_1 + (70 << Exp_shift);
993 else if (!oldinexact)
996 #ifdef Avoid_Underflow
998 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
1000 dval(rv) *= dval(rv0);
1002 /* try to avoid the bug of testing an 8087 register value */
1004 if (!(word0(rv) & Exp_mask))
1006 if (word0(rv) == 0 && word1(rv) == 0)
1011 #endif /* Avoid_Underflow */
1013 if (inexact && !(word0(rv) & Exp_mask)) {
1014 /* set underflow bit */
1016 dval(rv0) *= dval(rv0);
1028 return sign ? -dval(rv) : dval(rv);