Deleted Added
sdiff udiff text old ( 187808 ) new ( 219557 )
full compact
strtodg.c (187808) strtodg.c (219557)
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998-2001 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and

--- 158 unchanged lines hidden (view full) ---

167 x[-1] >>= ULbits - n;
168 return b;
169 }
170
171 static int
172rvOK
173#ifdef KR_headers
174 (d, fpi, exp, bits, exact, rd, irv)
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998-2001 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and

--- 158 unchanged lines hidden (view full) ---

167 x[-1] >>= ULbits - n;
168 return b;
169 }
170
171 static int
172rvOK
173#ifdef KR_headers
174 (d, fpi, exp, bits, exact, rd, irv)
175 double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
175 U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
176#else
176#else
177 (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
177 (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
178#endif
179{
180 Bigint *b;
181 ULong carry, inex, lostbits;
182 int bdif, e, j, k, k1, nb, rv;
183
184 carry = rv = 0;
178#endif
179{
180 Bigint *b;
181 ULong carry, inex, lostbits;
182 int bdif, e, j, k, k1, nb, rv;
183
184 carry = rv = 0;
185 b = d2b(d, &e, &bdif);
185 b = d2b(dval(d), &e, &bdif);
186 bdif -= nb = fpi->nbits;
187 e += bdif;
188 if (bdif <= 0) {
189 if (exact)
190 goto trunc;
191 goto ret;
192 }
193 if (P == nb) {

--- 92 unchanged lines hidden (view full) ---

286 rv = 1;
287 ret:
288 Bfree(b);
289 return rv;
290 }
291
292 static int
293#ifdef KR_headers
186 bdif -= nb = fpi->nbits;
187 e += bdif;
188 if (bdif <= 0) {
189 if (exact)
190 goto trunc;
191 goto ret;
192 }
193 if (P == nb) {

--- 92 unchanged lines hidden (view full) ---

286 rv = 1;
287 ret:
288 Bfree(b);
289 return rv;
290 }
291
292 static int
293#ifdef KR_headers
294mantbits(d) double d;
294mantbits(d) U *d;
295#else
295#else
296mantbits(double d)
296mantbits(U *d)
297#endif
298{
299 ULong L;
300#ifdef VAX
301 L = word1(d) << 16 | word1(d) >> 16;
302 if (L)
303#else
304 if ( (L = word1(d)) !=0)

--- 17 unchanged lines hidden (view full) ---

322#endif
323{
324 int abe, abits, asub;
325 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
326 int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
327 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
328 int sudden_underflow;
329 CONST char *s, *s0, *s1;
297#endif
298{
299 ULong L;
300#ifdef VAX
301 L = word1(d) << 16 | word1(d) >> 16;
302 if (L)
303#else
304 if ( (L = word1(d)) !=0)

--- 17 unchanged lines hidden (view full) ---

322#endif
323{
324 int abe, abits, asub;
325 int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
326 int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
327 int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
328 int sudden_underflow;
329 CONST char *s, *s0, *s1;
330 double adj, adj0, rv, tol;
330 double adj0, tol;
331 Long L;
331 Long L;
332 U adj, rv;
332 ULong *b, *be, y, z;
333 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
334#ifdef USE_LOCALE /*{{*/
335#ifdef NO_LOCALE_CACHE
336 char *decimalpoint = localeconv()->decimal_point;
337 int dplen = strlen(decimalpoint);
338#else
339 char *decimalpoint;
340 static char *decimalpoint_cache;
341 static int dplen;
342 if (!(s0 = decimalpoint_cache)) {
343 s0 = localeconv()->decimal_point;
333 ULong *b, *be, y, z;
334 Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
335#ifdef USE_LOCALE /*{{*/
336#ifdef NO_LOCALE_CACHE
337 char *decimalpoint = localeconv()->decimal_point;
338 int dplen = strlen(decimalpoint);
339#else
340 char *decimalpoint;
341 static char *decimalpoint_cache;
342 static int dplen;
343 if (!(s0 = decimalpoint_cache)) {
344 s0 = localeconv()->decimal_point;
344 if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
345 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
345 strcpy(decimalpoint_cache, s0);
346 s0 = decimalpoint_cache;
347 }
348 dplen = strlen(s0);
349 }
350 decimalpoint = (char*)s0;
351#endif /*NO_LOCALE_CACHE*/
352#else /*USE_LOCALE}{*/
353#define dplen 1
354#endif /*USE_LOCALE}}*/
355
356 irv = STRTOG_Zero;
357 denorm = sign = nz0 = nz = 0;
346 strcpy(decimalpoint_cache, s0);
347 s0 = decimalpoint_cache;
348 }
349 dplen = strlen(s0);
350 }
351 decimalpoint = (char*)s0;
352#endif /*NO_LOCALE_CACHE*/
353#else /*USE_LOCALE}{*/
354#define dplen 1
355#endif /*USE_LOCALE}}*/
356
357 irv = STRTOG_Zero;
358 denorm = sign = nz0 = nz = 0;
358 dval(rv) = 0.;
359 dval(&rv) = 0.;
359 rvb = 0;
360 nbits = fpi->nbits;
361 for(s = s00;;s++) switch(*s) {
362 case '-':
363 sign = 1;
364 /* no break */
365 case '+':
366 if (*++s)

--- 175 unchanged lines hidden (view full) ---

542 /* Now we have nd0 digits, starting at s0, followed by a
543 * decimal point, followed by nd-nd0 digits. The number we're
544 * after is the integer represented by those digits times
545 * 10**e */
546
547 if (!nd0)
548 nd0 = nd;
549 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
360 rvb = 0;
361 nbits = fpi->nbits;
362 for(s = s00;;s++) switch(*s) {
363 case '-':
364 sign = 1;
365 /* no break */
366 case '+':
367 if (*++s)

--- 175 unchanged lines hidden (view full) ---

543 /* Now we have nd0 digits, starting at s0, followed by a
544 * decimal point, followed by nd-nd0 digits. The number we're
545 * after is the integer represented by those digits times
546 * 10**e */
547
548 if (!nd0)
549 nd0 = nd;
550 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
550 dval(rv) = y;
551 dval(&rv) = y;
551 if (k > 9)
552 if (k > 9)
552 dval(rv) = tens[k - 9] * dval(rv) + z;
553 dval(&rv) = tens[k - 9] * dval(&rv) + z;
553 bd0 = 0;
554 if (nbits <= P && nd <= DBL_DIG) {
555 if (!e) {
554 bd0 = 0;
555 if (nbits <= P && nd <= DBL_DIG) {
556 if (!e) {
556 if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
557 if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
557 goto ret;
558 }
559 else if (e > 0) {
560 if (e <= Ten_pmax) {
561#ifdef VAX
562 goto vax_ovfl_check;
563#else
558 goto ret;
559 }
560 else if (e > 0) {
561 if (e <= Ten_pmax) {
562#ifdef VAX
563 goto vax_ovfl_check;
564#else
564 i = fivesbits[e] + mantbits(dval(rv)) <= P;
565 /* rv = */ rounded_product(dval(rv), tens[e]);
566 if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
565 i = fivesbits[e] + mantbits(&rv) <= P;
566 /* rv = */ rounded_product(dval(&rv), tens[e]);
567 if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
567 goto ret;
568 e1 -= e;
569 goto rv_notOK;
570#endif
571 }
572 i = DBL_DIG - nd;
573 if (e <= Ten_pmax + i) {
574 /* A fancier test would sometimes let us do
575 * this for larger i values.
576 */
577 e2 = e - i;
578 e1 -= i;
568 goto ret;
569 e1 -= e;
570 goto rv_notOK;
571#endif
572 }
573 i = DBL_DIG - nd;
574 if (e <= Ten_pmax + i) {
575 /* A fancier test would sometimes let us do
576 * this for larger i values.
577 */
578 e2 = e - i;
579 e1 -= i;
579 dval(rv) *= tens[i];
580 dval(&rv) *= tens[i];
580#ifdef VAX
581 /* VAX exponent range is so narrow we must
582 * worry about overflow here...
583 */
584 vax_ovfl_check:
581#ifdef VAX
582 /* VAX exponent range is so narrow we must
583 * worry about overflow here...
584 */
585 vax_ovfl_check:
585 dval(adj) = dval(rv);
586 word0(adj) -= P*Exp_msk1;
587 /* adj = */ rounded_product(dval(adj), tens[e2]);
588 if ((word0(adj) & Exp_mask)
586 dval(&adj) = dval(&rv);
587 word0(&adj) -= P*Exp_msk1;
588 /* adj = */ rounded_product(dval(&adj), tens[e2]);
589 if ((word0(&adj) & Exp_mask)
589 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
590 goto rv_notOK;
590 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
591 goto rv_notOK;
591 word0(adj) += P*Exp_msk1;
592 dval(rv) = dval(adj);
592 word0(&adj) += P*Exp_msk1;
593 dval(&rv) = dval(&adj);
593#else
594#else
594 /* rv = */ rounded_product(dval(rv), tens[e2]);
595 /* rv = */ rounded_product(dval(&rv), tens[e2]);
595#endif
596#endif
596 if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
597 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
597 goto ret;
598 e1 -= e2;
599 }
600 }
601#ifndef Inaccurate_Divide
602 else if (e >= -Ten_pmax) {
598 goto ret;
599 e1 -= e2;
600 }
601 }
602#ifndef Inaccurate_Divide
603 else if (e >= -Ten_pmax) {
603 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
604 if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
604 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
605 if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
605 goto ret;
606 e1 -= e;
607 }
608#endif
609 }
610 rv_notOK:
611 e1 += nd - k;
612
613 /* Get starting approximation = rv * 10**e1 */
614
615 e2 = 0;
616 if (e1 > 0) {
617 if ( (i = e1 & 15) !=0)
606 goto ret;
607 e1 -= e;
608 }
609#endif
610 }
611 rv_notOK:
612 e1 += nd - k;
613
614 /* Get starting approximation = rv * 10**e1 */
615
616 e2 = 0;
617 if (e1 > 0) {
618 if ( (i = e1 & 15) !=0)
618 dval(rv) *= tens[i];
619 dval(&rv) *= tens[i];
619 if (e1 &= ~15) {
620 e1 >>= 4;
620 if (e1 &= ~15) {
621 e1 >>= 4;
621 while(e1 >= (1 << n_bigtens-1)) {
622 e2 += ((word0(rv) & Exp_mask)
622 while(e1 >= (1 << (n_bigtens-1))) {
623 e2 += ((word0(&rv) & Exp_mask)
623 >> Exp_shift1) - Bias;
624 >> Exp_shift1) - Bias;
624 word0(rv) &= ~Exp_mask;
625 word0(rv) |= Bias << Exp_shift1;
626 dval(rv) *= bigtens[n_bigtens-1];
627 e1 -= 1 << n_bigtens-1;
625 word0(&rv) &= ~Exp_mask;
626 word0(&rv) |= Bias << Exp_shift1;
627 dval(&rv) *= bigtens[n_bigtens-1];
628 e1 -= 1 << (n_bigtens-1);
628 }
629 }
629 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
630 word0(rv) &= ~Exp_mask;
631 word0(rv) |= Bias << Exp_shift1;
630 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
631 word0(&rv) &= ~Exp_mask;
632 word0(&rv) |= Bias << Exp_shift1;
632 for(j = 0; e1 > 0; j++, e1 >>= 1)
633 if (e1 & 1)
633 for(j = 0; e1 > 0; j++, e1 >>= 1)
634 if (e1 & 1)
634 dval(rv) *= bigtens[j];
635 dval(&rv) *= bigtens[j];
635 }
636 }
637 else if (e1 < 0) {
638 e1 = -e1;
639 if ( (i = e1 & 15) !=0)
636 }
637 }
638 else if (e1 < 0) {
639 e1 = -e1;
640 if ( (i = e1 & 15) !=0)
640 dval(rv) /= tens[i];
641 dval(&rv) /= tens[i];
641 if (e1 &= ~15) {
642 e1 >>= 4;
642 if (e1 &= ~15) {
643 e1 >>= 4;
643 while(e1 >= (1 << n_bigtens-1)) {
644 e2 += ((word0(rv) & Exp_mask)
644 while(e1 >= (1 << (n_bigtens-1))) {
645 e2 += ((word0(&rv) & Exp_mask)
645 >> Exp_shift1) - Bias;
646 >> Exp_shift1) - Bias;
646 word0(rv) &= ~Exp_mask;
647 word0(rv) |= Bias << Exp_shift1;
648 dval(rv) *= tinytens[n_bigtens-1];
649 e1 -= 1 << n_bigtens-1;
647 word0(&rv) &= ~Exp_mask;
648 word0(&rv) |= Bias << Exp_shift1;
649 dval(&rv) *= tinytens[n_bigtens-1];
650 e1 -= 1 << (n_bigtens-1);
650 }
651 }
651 e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
652 word0(rv) &= ~Exp_mask;
653 word0(rv) |= Bias << Exp_shift1;
652 e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
653 word0(&rv) &= ~Exp_mask;
654 word0(&rv) |= Bias << Exp_shift1;
654 for(j = 0; e1 > 0; j++, e1 >>= 1)
655 if (e1 & 1)
655 for(j = 0; e1 > 0; j++, e1 >>= 1)
656 if (e1 & 1)
656 dval(rv) *= tinytens[j];
657 dval(&rv) *= tinytens[j];
657 }
658 }
659#ifdef IBM
660 /* e2 is a correction to the (base 2) exponent of the return
661 * value, reflecting adjustments above to avoid overflow in the
662 * native arithmetic. For native IBM (base 16) arithmetic, we
663 * must multiply e2 by 4 to change from base 16 to 2.
664 */
665 e2 <<= 2;
666#endif
658 }
659 }
660#ifdef IBM
661 /* e2 is a correction to the (base 2) exponent of the return
662 * value, reflecting adjustments above to avoid overflow in the
663 * native arithmetic. For native IBM (base 16) arithmetic, we
664 * must multiply e2 by 4 to change from base 16 to 2.
665 */
666 e2 <<= 2;
667#endif
667 rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */
668 rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
668 rve += e2;
669 if ((j = rvbits - nbits) > 0) {
670 rshift(rvb, j);
671 rvbits = nbits;
672 rve += j;
673 }
674 bb0 = 0; /* trailing zero bits in rvb */
675 e2 = rve + rvbits - nbits;

--- 161 unchanged lines hidden (view full) ---

837 break;
838 }
839 rve -= nbits;
840 rvb = set_ones(rvb, rvbits = nbits);
841 break;
842 }
843 else
844 irv = STRTOG_Normal | STRTOG_Inexhi;
669 rve += e2;
670 if ((j = rvbits - nbits) > 0) {
671 rshift(rvb, j);
672 rvbits = nbits;
673 rve += j;
674 }
675 bb0 = 0; /* trailing zero bits in rvb */
676 e2 = rve + rvbits - nbits;

--- 161 unchanged lines hidden (view full) ---

838 break;
839 }
840 rve -= nbits;
841 rvb = set_ones(rvb, rvbits = nbits);
842 break;
843 }
844 else
845 irv = STRTOG_Normal | STRTOG_Inexhi;
845 if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
846 if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
846 break;
847 if (dsign) {
848 rvb = increment(rvb);
849 j = kmask & (ULbits - (rvbits & kmask));
850 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
851 rvbits++;
852 irv = STRTOG_Normal | STRTOG_Inexhi;
853 }
854 else {
855 if (bbbits == 1)
856 goto undfl;
857 decrement(rvb);
858 irv = STRTOG_Normal | STRTOG_Inexlo;
859 }
860 break;
861 }
847 break;
848 if (dsign) {
849 rvb = increment(rvb);
850 j = kmask & (ULbits - (rvbits & kmask));
851 if (hi0bits(rvb->x[rvb->wds - 1]) != j)
852 rvbits++;
853 irv = STRTOG_Normal | STRTOG_Inexhi;
854 }
855 else {
856 if (bbbits == 1)
857 goto undfl;
858 decrement(rvb);
859 irv = STRTOG_Normal | STRTOG_Inexlo;
860 }
861 break;
862 }
862 if ((dval(adj) = ratio(delta, bs)) <= 2.) {
863 if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
863 adj1:
864 inex = STRTOG_Inexlo;
865 if (dsign) {
866 asub = 0;
867 inex = STRTOG_Inexhi;
868 }
869 else if (denorm && bbbits <= 1) {
870 undfl:
871 rvb->wds = 0;
872 rve = emin;
873 irv = STRTOG_Underflow | STRTOG_Inexlo;
874 break;
875 }
864 adj1:
865 inex = STRTOG_Inexlo;
866 if (dsign) {
867 asub = 0;
868 inex = STRTOG_Inexhi;
869 }
870 else if (denorm && bbbits <= 1) {
871 undfl:
872 rvb->wds = 0;
873 rve = emin;
874 irv = STRTOG_Underflow | STRTOG_Inexlo;
875 break;
876 }
876 adj0 = dval(adj) = 1.;
877 adj0 = dval(&adj) = 1.;
877 }
878 else {
878 }
879 else {
879 adj0 = dval(adj) *= 0.5;
880 adj0 = dval(&adj) *= 0.5;
880 if (dsign) {
881 asub = 0;
882 inex = STRTOG_Inexlo;
883 }
881 if (dsign) {
882 asub = 0;
883 inex = STRTOG_Inexlo;
884 }
884 if (dval(adj) < 2147483647.) {
885 if (dval(&adj) < 2147483647.) {
885 L = adj0;
886 adj0 -= L;
887 switch(rd) {
888 case 0:
889 if (adj0 >= .5)
890 goto inc_L;
891 break;
892 case 1:
893 if (asub && adj0 > 0.)
894 goto inc_L;
895 break;
896 case 2:
897 if (!asub && adj0 > 0.) {
898 inc_L:
899 L++;
900 inex = STRTOG_Inexact - inex;
901 }
902 }
886 L = adj0;
887 adj0 -= L;
888 switch(rd) {
889 case 0:
890 if (adj0 >= .5)
891 goto inc_L;
892 break;
893 case 1:
894 if (asub && adj0 > 0.)
895 goto inc_L;
896 break;
897 case 2:
898 if (!asub && adj0 > 0.) {
899 inc_L:
900 L++;
901 inex = STRTOG_Inexact - inex;
902 }
903 }
903 dval(adj) = L;
904 dval(&adj) = L;
904 }
905 }
906 y = rve + rvbits;
907
905 }
906 }
907 y = rve + rvbits;
908
908 /* adj *= ulp(dval(rv)); */
909 /* adj *= ulp(dval(&rv)); */
909 /* if (asub) rv -= adj; else rv += adj; */
910
911 if (!denorm && rvbits < nbits) {
912 rvb = lshift(rvb, j = nbits - rvbits);
913 rve -= j;
914 rvbits = nbits;
915 }
910 /* if (asub) rv -= adj; else rv += adj; */
911
912 if (!denorm && rvbits < nbits) {
913 rvb = lshift(rvb, j = nbits - rvbits);
914 rve -= j;
915 rvbits = nbits;
916 }
916 ab = d2b(dval(adj), &abe, &abits);
917 ab = d2b(dval(&adj), &abe, &abits);
917 if (abe < 0)
918 rshift(ab, -abe);
919 else if (abe > 0)
920 ab = lshift(ab, abe);
921 rvb0 = rvb;
922 if (asub) {
923 /* rv -= adj; */
924 j = hi0bits(rvb->x[rvb->wds-1]);

--- 37 unchanged lines hidden (view full) ---

962 Bfree(ab);
963 Bfree(rvb0);
964 if (finished)
965 break;
966
967 z = rve + rvbits;
968 if (y == z && L) {
969 /* Can we stop now? */
918 if (abe < 0)
919 rshift(ab, -abe);
920 else if (abe > 0)
921 ab = lshift(ab, abe);
922 rvb0 = rvb;
923 if (asub) {
924 /* rv -= adj; */
925 j = hi0bits(rvb->x[rvb->wds-1]);

--- 37 unchanged lines hidden (view full) ---

963 Bfree(ab);
964 Bfree(rvb0);
965 if (finished)
966 break;
967
968 z = rve + rvbits;
969 if (y == z && L) {
970 /* Can we stop now? */
970 tol = dval(adj) * 5e-16; /* > max rel error */
971 dval(adj) = adj0 - .5;
972 if (dval(adj) < -tol) {
971 tol = dval(&adj) * 5e-16; /* > max rel error */
972 dval(&adj) = adj0 - .5;
973 if (dval(&adj) < -tol) {
973 if (adj0 > tol) {
974 irv |= inex;
975 break;
976 }
977 }
974 if (adj0 > tol) {
975 irv |= inex;
976 break;
977 }
978 }
978 else if (dval(adj) > tol && adj0 < 1. - tol) {
979 else if (dval(&adj) > tol && adj0 < 1. - tol) {
979 irv |= inex;
980 break;
981 }
982 }
983 bb0 = denorm ? 0 : trailz(rvb);
984 Bfree(bb);
985 Bfree(bd);
986 Bfree(bs);

--- 78 unchanged lines hidden --- 		980 irv |= inex;
981 break;
982 }
983 }
984 bb0 = denorm ? 0 : trailz(rvb);
985 Bfree(bb);
986 Bfree(bd);
987 Bfree(bs);

--- 78 unchanged lines hidden ---