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
42 #define Avoid_Underflow
44 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
45 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
46 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
47 9007199254740992.e-256
52 #ifdef Honor_FLT_ROUNDS
53 #define Rounding rounding
54 #undef Check_FLT_ROUNDS
55 #define Check_FLT_ROUNDS
57 #define Rounding Flt_Rounds
63 (s00, se) CONST char *s00; char **se;
65 (CONST char *s00, char **se)
68 #ifdef Avoid_Underflow
71 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
72 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
73 CONST char *s, *s0, *s1;
74 double aadj, aadj1, adj, rv, rv0;
77 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
79 int inexact, oldinexact;
81 #ifdef Honor_FLT_ROUNDS
87 for(s = s00;;s++) switch(*s) {
111 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
117 switch((i = gethex(&s, &fpi, &exp, &bb, sign)) & STRTOG_Retmask) {
118 case STRTOG_NoNumber:
125 copybits(bits, fpi.nbits, bb);
128 ULtod(((U*)&rv)->L, bits, exp, i);
141 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
148 if (c == *localeconv()->decimal_point)
155 for(; c == '0'; c = *++s)
157 if (c > '0' && c <= '9') {
165 for(; c >= '0' && c <= '9'; c = *++s) {
170 for(i = 1; i < nz; i++)
173 else if (nd <= DBL_DIG + 1)
177 else if (nd <= DBL_DIG + 1)
185 if (c == 'e' || c == 'E') {
186 if (!nd && !nz && !nz0) {
197 if (c >= '0' && c <= '9') {
200 if (c > '0' && c <= '9') {
203 while((c = *++s) >= '0' && c <= '9')
205 if (s - s1 > 8 || L > 19999)
206 /* Avoid confusion from exponents
207 * so large that e might overflow.
209 e = 19999; /* safe for 16 bit ints */
224 /* Check for Nan and Infinity */
226 static FPI fpinan = /* only 52 explicit bits */
227 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
231 if (match(&s,"nf")) {
233 if (!match(&s,"inity"))
235 word0(rv) = 0x7ff00000;
242 if (match(&s, "an")) {
245 && hexnan(&s, &fpinan, bits)
247 word0(rv) = 0x7ff00000 | bits[1];
251 word0(rv) = NAN_WORD0;
252 word1(rv) = NAN_WORD1;
258 #endif /* INFNAN_CHECK */
267 /* Now we have nd0 digits, starting at s0, followed by a
268 * decimal point, followed by nd-nd0 digits. The number we're
269 * after is the integer represented by those digits times
274 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
279 oldinexact = get_inexact();
281 dval(rv) = tens[k - 9] * dval(rv) + z;
286 #ifndef Honor_FLT_ROUNDS
298 #ifdef Honor_FLT_ROUNDS
299 /* round correctly FLT_ROUNDS = 2 or 3 */
305 /* rv = */ rounded_product(dval(rv), tens[e]);
310 if (e <= Ten_pmax + i) {
311 /* A fancier test would sometimes let us do
312 * this for larger i values.
314 #ifdef Honor_FLT_ROUNDS
315 /* round correctly FLT_ROUNDS = 2 or 3 */
324 /* VAX exponent range is so narrow we must
325 * worry about overflow here...
328 word0(rv) -= P*Exp_msk1;
329 /* rv = */ rounded_product(dval(rv), tens[e]);
330 if ((word0(rv) & Exp_mask)
331 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
333 word0(rv) += P*Exp_msk1;
335 /* rv = */ rounded_product(dval(rv), tens[e]);
340 #ifndef Inaccurate_Divide
341 else if (e >= -Ten_pmax) {
342 #ifdef Honor_FLT_ROUNDS
343 /* round correctly FLT_ROUNDS = 2 or 3 */
349 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
360 oldinexact = get_inexact();
362 #ifdef Avoid_Underflow
365 #ifdef Honor_FLT_ROUNDS
366 if ((rounding = Flt_Rounds) >= 2) {
368 rounding = rounding == 2 ? 0 : 2;
374 #endif /*IEEE_Arith*/
376 /* Get starting approximation = rv * 10**e1 */
379 if ( (i = e1 & 15) !=0)
382 if (e1 > DBL_MAX_10_EXP) {
387 /* Can't trust HUGE_VAL */
389 #ifdef Honor_FLT_ROUNDS
391 case 0: /* toward 0 */
392 case 3: /* toward -infinity */
397 word0(rv) = Exp_mask;
400 #else /*Honor_FLT_ROUNDS*/
401 word0(rv) = Exp_mask;
403 #endif /*Honor_FLT_ROUNDS*/
405 /* set overflow bit */
407 dval(rv0) *= dval(rv0);
412 #endif /*IEEE_Arith*/
418 for(j = 0; e1 > 1; j++, e1 >>= 1)
420 dval(rv) *= bigtens[j];
421 /* The last multiplication could overflow. */
422 word0(rv) -= P*Exp_msk1;
423 dval(rv) *= bigtens[j];
424 if ((z = word0(rv) & Exp_mask)
425 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
427 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
428 /* set to largest number */
429 /* (Can't trust DBL_MAX) */
434 word0(rv) += P*Exp_msk1;
439 if ( (i = e1 & 15) !=0)
442 if (e1 >= 1 << n_bigtens)
444 #ifdef Avoid_Underflow
447 for(j = 0; e1 > 0; j++, e1 >>= 1)
449 dval(rv) *= tinytens[j];
450 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
451 >> Exp_shift)) > 0) {
452 /* scaled rv is denormal; zap j low bits */
456 word0(rv) = (P+2)*Exp_msk1;
458 word0(rv) &= 0xffffffff << j-32;
461 word1(rv) &= 0xffffffff << j;
464 for(j = 0; e1 > 1; j++, e1 >>= 1)
466 dval(rv) *= tinytens[j];
467 /* The last multiplication could underflow. */
468 dval(rv0) = dval(rv);
469 dval(rv) *= tinytens[j];
471 dval(rv) = 2.*dval(rv0);
472 dval(rv) *= tinytens[j];
484 #ifndef Avoid_Underflow
487 /* The refinement below will clean
488 * this approximation up.
495 /* Now the hard part -- adjusting rv to the correct value.*/
497 /* Put digits into bd: true value = bd * 10^e */
499 bd0 = s2b(s0, nd0, nd, y);
504 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
520 #ifdef Honor_FLT_ROUNDS
524 #ifdef Avoid_Underflow
526 i = j + bbbits - 1; /* logb(rv) */
527 if (i < Emin) /* denormal */
531 #else /*Avoid_Underflow*/
532 #ifdef Sudden_Underflow
534 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
538 #else /*Sudden_Underflow*/
540 i = j + bbbits - 1; /* logb(rv) */
541 if (i < Emin) /* denormal */
545 #endif /*Sudden_Underflow*/
546 #endif /*Avoid_Underflow*/
549 #ifdef Avoid_Underflow
552 i = bb2 < bd2 ? bb2 : bd2;
561 bs = pow5mult(bs, bb5);
567 bb = lshift(bb, bb2);
569 bd = pow5mult(bd, bd5);
571 bd = lshift(bd, bd2);
573 bs = lshift(bs, bs2);
574 delta = diff(bb, bd);
578 #ifdef Honor_FLT_ROUNDS
581 /* Error is less than an ulp */
582 if (!delta->x[0] && delta->wds <= 1) {
598 && !(word0(rv) & Frac_mask)) {
599 y = word0(rv) & Exp_mask;
600 #ifdef Avoid_Underflow
601 if (!scale || y > 2*P*Exp_msk1)
606 delta = lshift(delta,Log2P);
607 if (cmp(delta, bs) <= 0)
612 #ifdef Avoid_Underflow
613 if (scale && (y = word0(rv) & Exp_mask)
615 word0(adj) += (2*P+1)*Exp_msk1 - y;
617 #ifdef Sudden_Underflow
618 if ((word0(rv) & Exp_mask) <=
620 word0(rv) += P*Exp_msk1;
621 dval(rv) += adj*ulp(dval(rv));
622 word0(rv) -= P*Exp_msk1;
625 #endif /*Sudden_Underflow*/
626 #endif /*Avoid_Underflow*/
627 dval(rv) += adj*ulp(dval(rv));
631 adj = ratio(delta, bs);
634 if (adj <= 0x7ffffffe) {
635 /* adj = rounding ? ceil(adj) : floor(adj); */
638 if (!((rounding>>1) ^ dsign))
643 #ifdef Avoid_Underflow
644 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
645 word0(adj) += (2*P+1)*Exp_msk1 - y;
647 #ifdef Sudden_Underflow
648 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
649 word0(rv) += P*Exp_msk1;
650 adj *= ulp(dval(rv));
655 word0(rv) -= P*Exp_msk1;
658 #endif /*Sudden_Underflow*/
659 #endif /*Avoid_Underflow*/
660 adj *= ulp(dval(rv));
667 #endif /*Honor_FLT_ROUNDS*/
670 /* Error is less than half an ulp -- check for
671 * special case of mantissa a power of two.
673 if (dsign || word1(rv) || word0(rv) & Bndry_mask
675 #ifdef Avoid_Underflow
676 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
678 || (word0(rv) & Exp_mask) <= Exp_msk1
683 if (!delta->x[0] && delta->wds <= 1)
688 if (!delta->x[0] && delta->wds <= 1) {
695 delta = lshift(delta,Log2P);
696 if (cmp(delta, bs) > 0)
701 /* exactly half-way between */
703 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
705 #ifdef Avoid_Underflow
706 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
707 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
710 /*boundary case -- increment exponent*/
711 word0(rv) = (word0(rv) & Exp_mask)
718 #ifdef Avoid_Underflow
724 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
726 /* boundary case -- decrement exponent */
727 #ifdef Sudden_Underflow /*{{*/
728 L = word0(rv) & Exp_mask;
732 #ifdef Avoid_Underflow
733 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
736 #endif /*Avoid_Underflow*/
740 #else /*Sudden_Underflow}{*/
741 #ifdef Avoid_Underflow
743 L = word0(rv) & Exp_mask;
744 if (L <= (2*P+1)*Exp_msk1) {
745 if (L > (P+2)*Exp_msk1)
749 /* rv = smallest denormal */
753 #endif /*Avoid_Underflow*/
754 L = (word0(rv) & Exp_mask) - Exp_msk1;
755 #endif /*Sudden_Underflow}*/
756 word0(rv) = L | Bndry_mask1;
757 word1(rv) = 0xffffffff;
765 if (!(word1(rv) & LSB))
769 dval(rv) += ulp(dval(rv));
772 dval(rv) -= ulp(dval(rv));
773 #ifndef Sudden_Underflow
778 #ifdef Avoid_Underflow
784 if ((aadj = ratio(delta, bs)) <= 2.) {
787 else if (word1(rv) || word0(rv) & Bndry_mask) {
788 #ifndef Sudden_Underflow
789 if (word1(rv) == Tiny1 && !word0(rv))
796 /* special case -- power of FLT_RADIX to be */
797 /* rounded down... */
799 if (aadj < 2./FLT_RADIX)
808 aadj1 = dsign ? aadj : -aadj;
809 #ifdef Check_FLT_ROUNDS
811 case 2: /* towards +infinity */
814 case 0: /* towards 0 */
815 case 3: /* towards -infinity */
821 #endif /*Check_FLT_ROUNDS*/
823 y = word0(rv) & Exp_mask;
825 /* Check for overflow */
827 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
828 dval(rv0) = dval(rv);
829 word0(rv) -= P*Exp_msk1;
830 adj = aadj1 * ulp(dval(rv));
832 if ((word0(rv) & Exp_mask) >=
833 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
834 if (word0(rv0) == Big0 && word1(rv0) == Big1)
841 word0(rv) += P*Exp_msk1;
844 #ifdef Avoid_Underflow
845 if (scale && y <= 2*P*Exp_msk1) {
846 if (aadj <= 0x7fffffff) {
850 aadj1 = dsign ? aadj : -aadj;
852 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
854 adj = aadj1 * ulp(dval(rv));
857 #ifdef Sudden_Underflow
858 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
859 dval(rv0) = dval(rv);
860 word0(rv) += P*Exp_msk1;
861 adj = aadj1 * ulp(dval(rv));
864 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
866 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
869 if (word0(rv0) == Tiny0
870 && word1(rv0) == Tiny1)
877 word0(rv) -= P*Exp_msk1;
880 adj = aadj1 * ulp(dval(rv));
883 #else /*Sudden_Underflow*/
884 /* Compute adj so that the IEEE rounding rules will
885 * correctly round rv + adj in some half-way cases.
886 * If rv * ulp(rv) is denormalized (i.e.,
887 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
888 * trouble from bits lost to denormalization;
889 * example: 1.2e-307 .
891 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
892 aadj1 = (double)(int)(aadj + 0.5);
896 adj = aadj1 * ulp(dval(rv));
898 #endif /*Sudden_Underflow*/
899 #endif /*Avoid_Underflow*/
901 z = word0(rv) & Exp_mask;
903 #ifdef Avoid_Underflow
907 /* Can we stop now? */
910 /* The tolerances below are conservative. */
911 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
912 if (aadj < .4999999 || aadj > .5000001)
915 else if (aadj < .4999999/FLT_RADIX)
928 word0(rv0) = Exp_1 + (70 << Exp_shift);
933 else if (!oldinexact)
936 #ifdef Avoid_Underflow
938 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
940 dval(rv) *= dval(rv0);
942 /* try to avoid the bug of testing an 8087 register value */
943 if (word0(rv) == 0 && word1(rv) == 0)
947 #endif /* Avoid_Underflow */
949 if (inexact && !(word0(rv) & Exp_mask)) {
950 /* set underflow bit */
952 dval(rv0) *= dval(rv0);
964 return sign ? -dval(rv) : dval(rv);