Import gdtoa-20101105.
authorPeter Avalos <pavalos@dragonflybsd.org>
Mon, 31 Jan 2011 05:00:55 +0000 (19:00 -1000)
committerPeter Avalos <pavalos@dragonflybsd.org>
Mon, 31 Jan 2011 05:00:55 +0000 (19:00 -1000)
26 files changed:
contrib/gdtoa/README
contrib/gdtoa/dtoa.c
contrib/gdtoa/g__fmt.c
contrib/gdtoa/g_ddfmt.c
contrib/gdtoa/g_dfmt.c
contrib/gdtoa/gdtoa.c
contrib/gdtoa/gdtoaimp.h
contrib/gdtoa/gethex.c
contrib/gdtoa/hexnan.c
contrib/gdtoa/makefile
contrib/gdtoa/misc.c
contrib/gdtoa/smisc.c
contrib/gdtoa/strtoIg.c
contrib/gdtoa/strtod.c
contrib/gdtoa/strtodI.c
contrib/gdtoa/strtodg.c
contrib/gdtoa/strtof.c
contrib/gdtoa/strtopdd.c
contrib/gdtoa/strtopf.c
contrib/gdtoa/strtopx.c
contrib/gdtoa/strtopxL.c
contrib/gdtoa/strtordd.c
contrib/gdtoa/strtorf.c
contrib/gdtoa/strtorx.c
contrib/gdtoa/strtorxL.c
contrib/gdtoa/ulp.c

index ce8be55..0c034f1 100644 (file)
@@ -353,5 +353,12 @@ you also compile with -DNO_LOCALE_CACHE, the details about the
 current "decimal point" character string are cached and assumed not
 to change during the program's execution.
 
+On machines with a 64-bit long double and perhaps a 113-bit "quad"
+type, you can invoke "make Printf" to add Printf (and variants, such
+as Fprintf) to gdtoa.a.  These are analogs, declared in stdio1.h, of
+printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
+double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
+precision printing.
+
 Please send comments to        David M. Gay (dmg at acm dot org, with " at "
 changed at "@" and " dot " changed to ".").
index 48fdf5e..6e5de30 100644 (file)
@@ -75,10 +75,10 @@ THIS SOFTWARE.
  char *
 dtoa
 #ifdef KR_headers
-       (d, mode, ndigits, decpt, sign, rve)
-       double d; int mode, ndigits, *decpt, *sign; char **rve;
+       (d0, mode, ndigits, decpt, sign, rve)
+       double d0; int mode, ndigits, *decpt, *sign; char **rve;
 #else
-       (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
+       (double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
 #endif
 {
  /*    Arguments ndigits, decpt, sign are similar to those
@@ -124,7 +124,8 @@ dtoa
        ULong x;
 #endif
        Bigint *b, *b1, *delta, *mlo, *mhi, *S;
-       double d2, ds, eps;
+       U d, d2, eps;
+       double ds;
        char *s, *s0;
 #ifdef SET_INEXACT
        int inexact, oldinexact;
@@ -149,35 +150,35 @@ dtoa
                dtoa_result = 0;
                }
 #endif
-
-       if (word0(d) & Sign_bit) {
+       d.d = d0;
+       if (word0(&d) & Sign_bit) {
                /* set sign for everything, including 0's and NaNs */
                *sign = 1;
-               word0(d) &= ~Sign_bit;  /* clear sign bit */
+               word0(&d) &= ~Sign_bit; /* clear sign bit */
                }
        else
                *sign = 0;
 
 #if defined(IEEE_Arith) + defined(VAX)
 #ifdef IEEE_Arith
-       if ((word0(d) & Exp_mask) == Exp_mask)
+       if ((word0(&d) & Exp_mask) == Exp_mask)
 #else
-       if (word0(d)  == 0x8000)
+       if (word0(&d)  == 0x8000)
 #endif
                {
                /* Infinity or NaN */
                *decpt = 9999;
 #ifdef IEEE_Arith
-               if (!word1(d) && !(word0(d) & 0xfffff))
+               if (!word1(&d) && !(word0(&d) & 0xfffff))
                        return nrv_alloc("Infinity", rve, 8);
 #endif
                return nrv_alloc("NaN", rve, 3);
                }
 #endif
 #ifdef IBM
-       dval(d) += 0; /* normalize */
+       dval(&d) += 0; /* normalize */
 #endif
-       if (!dval(d)) {
+       if (!dval(&d)) {
                *decpt = 1;
                return nrv_alloc("0", rve, 1);
                }
@@ -196,26 +197,26 @@ dtoa
                }
 #endif
 
-       b = d2b(dval(d), &be, &bbits);
+       b = d2b(dval(&d), &be, &bbits);
 #ifdef Sudden_Underflow
-       i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
+       i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 #else
-       if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
+       if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
 #endif
-               dval(d2) = dval(d);
-               word0(d2) &= Frac_mask1;
-               word0(d2) |= Exp_11;
+               dval(&d2) = dval(&d);
+               word0(&d2) &= Frac_mask1;
+               word0(&d2) |= Exp_11;
 #ifdef IBM
-               if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
-                       dval(d2) /= 1 << j;
+               if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
+                       dval(&d2) /= 1 << j;
 #endif
 
                /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
                 * log10(x)      =  log(x) / log(10)
                 *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
-                * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
+                * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
                 *
-                * This suggests computing an approximation k to log10(d) by
+                * This suggests computing an approximation k to log10(&d) by
                 *
                 * k = (i - Bias)*0.301029995663981
                 *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@@ -244,21 +245,21 @@ dtoa
                /* d is denormalized */
 
                i = bbits + be + (Bias + (P-1) - 1);
-               x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
-                           : word1(d) << 32 - i;
-               dval(d2) = x;
-               word0(d2) -= 31*Exp_msk1; /* adjust exponent */
+               x = i > 32  ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
+                           : word1(&d) << (32 - i);
+               dval(&d2) = x;
+               word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
                i -= (Bias + (P-1) - 1) + 1;
                denorm = 1;
                }
 #endif
-       ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
+       ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
        k = (int)ds;
        if (ds < 0. && ds != k)
                k--;    /* want k = floor(ds) */
        k_check = 1;
        if (k >= 0 && k <= Ten_pmax) {
-               if (dval(d) < tens[k])
+               if (dval(&d) < tens[k])
                        k--;
                k_check = 0;
                }
@@ -297,10 +298,11 @@ dtoa
                try_quick = 0;
                }
        leftright = 1;
+       ilim = ilim1 = -1;      /* Values for cases 0 and 1; done here to */
+                               /* silence erroneous "gcc -Wall" warning. */
        switch(mode) {
                case 0:
                case 1:
-                       ilim = ilim1 = -1;
                        i = 18;
                        ndigits = 0;
                        break;
@@ -334,7 +336,7 @@ dtoa
                /* Try to get by with floating-point arithmetic. */
 
                i = 0;
-               dval(d2) = dval(d);
+               dval(&d2) = dval(&d);
                k0 = k;
                ilim0 = ilim;
                ieps = 2; /* conservative */
@@ -344,7 +346,7 @@ dtoa
                        if (j & Bletch) {
                                /* prevent overflows */
                                j &= Bletch - 1;
-                               dval(d) /= bigtens[n_bigtens-1];
+                               dval(&d) /= bigtens[n_bigtens-1];
                                ieps++;
                                }
                        for(; j; j >>= 1, i++)
@@ -352,32 +354,32 @@ dtoa
                                        ieps++;
                                        ds *= bigtens[i];
                                        }
-                       dval(d) /= ds;
+                       dval(&d) /= ds;
                        }
                else if (( j1 = -k )!=0) {
-                       dval(d) *= tens[j1 & 0xf];
+                       dval(&d) *= tens[j1 & 0xf];
                        for(j = j1 >> 4; j; j >>= 1, i++)
                                if (j & 1) {
                                        ieps++;
-                                       dval(d) *= bigtens[i];
+                                       dval(&d) *= bigtens[i];
                                        }
                        }
-               if (k_check && dval(d) < 1. && ilim > 0) {
+               if (k_check && dval(&d) < 1. && ilim > 0) {
                        if (ilim1 <= 0)
                                goto fast_failed;
                        ilim = ilim1;
                        k--;
-                       dval(d) *= 10.;
+                       dval(&d) *= 10.;
                        ieps++;
                        }
-               dval(eps) = ieps*dval(d) + 7.;
-               word0(eps) -= (P-1)*Exp_msk1;
+               dval(&eps) = ieps*dval(&d) + 7.;
+               word0(&eps) -= (P-1)*Exp_msk1;
                if (ilim == 0) {
                        S = mhi = 0;
-                       dval(d) -= 5.;
-                       if (dval(d) > dval(eps))
+                       dval(&d) -= 5.;
+                       if (dval(&d) > dval(&eps))
                                goto one_digit;
-                       if (dval(d) < -dval(eps))
+                       if (dval(&d) < -dval(&eps))
                                goto no_digits;
                        goto fast_failed;
                        }
@@ -386,34 +388,34 @@ dtoa
                        /* Use Steele & White method of only
                         * generating digits needed.
                         */
-                       dval(eps) = 0.5/tens[ilim-1] - dval(eps);
+                       dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
                        for(i = 0;;) {
-                               L = dval(d);
-                               dval(d) -= L;
+                               L = dval(&d);
+                               dval(&d) -= L;
                                *s++ = '0' + (int)L;
-                               if (dval(d) < dval(eps))
+                               if (dval(&d) < dval(&eps))
                                        goto ret1;
-                               if (1. - dval(d) < dval(eps))
+                               if (1. - dval(&d) < dval(&eps))
                                        goto bump_up;
                                if (++i >= ilim)
                                        break;
-                               dval(eps) *= 10.;
-                               dval(d) *= 10.;
+                               dval(&eps) *= 10.;
+                               dval(&d) *= 10.;
                                }
                        }
                else {
 #endif
                        /* Generate ilim digits, then fix them up. */
-                       dval(eps) *= tens[ilim-1];
-                       for(i = 1;; i++, dval(d) *= 10.) {
-                               L = (Long)(dval(d));
-                               if (!(dval(d) -= L))
+                       dval(&eps) *= tens[ilim-1];
+                       for(i = 1;; i++, dval(&d) *= 10.) {
+                               L = (Long)(dval(&d));
+                               if (!(dval(&d) -= L))
                                        ilim = i;
                                *s++ = '0' + (int)L;
                                if (i == ilim) {
-                                       if (dval(d) > 0.5 + dval(eps))
+                                       if (dval(&d) > 0.5 + dval(&eps))
                                                goto bump_up;
-                                       else if (dval(d) < 0.5 - dval(eps)) {
+                                       else if (dval(&d) < 0.5 - dval(&eps)) {
                                                while(*--s == '0');
                                                s++;
                                                goto ret1;
@@ -426,7 +428,7 @@ dtoa
 #endif
  fast_failed:
                s = s0;
-               dval(d) = dval(d2);
+               dval(&d) = dval(&d2);
                k = k0;
                ilim = ilim0;
                }
@@ -438,22 +440,22 @@ dtoa
                ds = tens[k];
                if (ndigits < 0 && ilim <= 0) {
                        S = mhi = 0;
-                       if (ilim < 0 || dval(d) <= 5*ds)
+                       if (ilim < 0 || dval(&d) <= 5*ds)
                                goto no_digits;
                        goto one_digit;
                        }
-               for(i = 1;; i++, dval(d) *= 10.) {
-                       L = (Long)(dval(d) / ds);
-                       dval(d) -= L*ds;
+               for(i = 1;; i++, dval(&d) *= 10.) {
+                       L = (Long)(dval(&d) / ds);
+                       dval(&d) -= L*ds;
 #ifdef Check_FLT_ROUNDS
                        /* If FLT_ROUNDS == 2, L will usually be high by 1 */
-                       if (dval(d) < 0) {
+                       if (dval(&d) < 0) {
                                L--;
-                               dval(d) += ds;
+                               dval(&d) += ds;
                                }
 #endif
                        *s++ = '0' + (int)L;
-                       if (!dval(d)) {
+                       if (!dval(&d)) {
 #ifdef SET_INEXACT
                                inexact = 0;
 #endif
@@ -467,8 +469,13 @@ dtoa
                                  case 2: goto bump_up;
                                  }
 #endif
-                               dval(d) += dval(d);
-                               if (dval(d) > ds || dval(d) == ds && L & 1) {
+                               dval(&d) += dval(&d);
+#ifdef ROUND_BIASED
+                               if (dval(&d) >= ds)
+#else
+                               if (dval(&d) > ds || (dval(&d) == ds && L & 1))
+#endif
+                                       {
  bump_up:
                                        while(*--s == '9')
                                                if (s == s0) {
@@ -533,9 +540,9 @@ dtoa
                        && Rounding == 1
 #endif
                                ) {
-               if (!word1(d) && !(word0(d) & Bndry_mask)
+               if (!word1(&d) && !(word0(&d) & Bndry_mask)
 #ifndef Sudden_Underflow
-                && word0(d) & (Exp_mask & ~Exp_msk1)
+                && word0(&d) & (Exp_mask & ~Exp_msk1)
 #endif
                                ) {
                        /* The special case */
@@ -621,7 +628,7 @@ dtoa
                        j1 = delta->sign ? 1 : cmp(b, delta);
                        Bfree(delta);
 #ifndef ROUND_BIASED
-                       if (j1 == 0 && mode != 1 && !(word1(d) & 1)
+                       if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
 #ifdef Honor_FLT_ROUNDS
                                && Rounding >= 1
 #endif
@@ -638,11 +645,11 @@ dtoa
                                goto ret;
                                }
 #endif
-                       if (j < 0 || j == 0 && mode != 1
+                       if (j < 0 || (j == 0 && mode != 1
 #ifndef ROUND_BIASED
-                                                       && !(word1(d) & 1)
+                                                       && !(word1(&d) & 1)
 #endif
-                                       ) {
+                                       )) {
                                if (!b->x[0] && b->wds <= 1) {
 #ifdef SET_INEXACT
                                        inexact = 0;
@@ -659,7 +666,11 @@ dtoa
                                if (j1 > 0) {
                                        b = lshift(b, 1);
                                        j1 = cmp(b, S);
-                                       if ((j1 > 0 || j1 == 0 && dig & 1)
+#ifdef ROUND_BIASED
+                                       if (j1 >= 0 /*)*/
+#else
+                                       if ((j1 > 0 || (j1 == 0 && dig & 1))
+#endif
                                        && dig++ == '9')
                                                goto round_9_up;
                                        }
@@ -719,7 +730,12 @@ dtoa
 #endif
        b = lshift(b, 1);
        j = cmp(b, S);
-       if (j > 0 || j == 0 && dig & 1) {
+#ifdef ROUND_BIASED
+       if (j >= 0)
+#else
+       if (j > 0 || (j == 0 && dig & 1))
+#endif
+               {
  roundoff:
                while(*--s == '9')
                        if (s == s0) {
@@ -730,7 +746,9 @@ dtoa
                ++*s++;
                }
        else {
+#ifdef Honor_FLT_ROUNDS
  trimzeros:
+#endif
                while(*--s == '0');
                s++;
                }
@@ -745,9 +763,9 @@ dtoa
 #ifdef SET_INEXACT
        if (inexact) {
                if (!oldinexact) {
-                       word0(d) = Exp_1 + (70 << Exp_shift);
-                       word1(d) = 0;
-                       dval(d) += 1.;
+                       word0(&d) = Exp_1 + (70 << Exp_shift);
+                       word1(&d) = 0;
+                       dval(&d) += 1.;
                        }
                }
        else if (!oldinexact)
index fbccb7d..6c4b3d2 100644 (file)
@@ -56,7 +56,7 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
        if (!(s0 = decimalpoint_cache)) {
                s0 = localeconv()->decimal_point;
                dlen = strlen(s0);
-               if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
+               if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
                        strcpy(decimalpoint_cache, s0);
                        s0 = decimalpoint_cache;
                        }
index b65d39d..1256655 100644 (file)
@@ -33,9 +33,9 @@ THIS SOFTWARE.
 
  char *
 #ifdef KR_headers
-g_ddfmt(buf, dd, ndig, bufsize) char *buf; double *dd; int ndig; size_t bufsize;
+g_ddfmt(buf, dd0, ndig, bufsize) char *buf; double *dd0; int ndig; size_t bufsize;
 #else
-g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
+g_ddfmt(char *buf, double *dd0, int ndig, size_t bufsize)
 #endif
 {
        FPI fpi;
@@ -43,7 +43,7 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
        ULong *L, bits0[4], *bits, *zx;
        int bx, by, decpt, ex, ey, i, j, mode;
        Bigint *x, *y, *z;
-       double ddx[2];
+       U *dd, ddx[2];
 #ifdef Honor_FLT_ROUNDS /*{{*/
        int Rounding;
 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
@@ -63,7 +63,8 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
        if (bufsize < 10 || bufsize < ndig + 8)
                return 0;
 
-       L = (ULong*)dd;
+       dd = (U*)dd0;
+       L = dd->L;
        if ((L[_0] & 0x7ff00000L) == 0x7ff00000L) {
                /* Infinity or NaN */
                if (L[_0] & 0xfffff || L[_1]) {
@@ -88,7 +89,7 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
                        goto nanret;
                goto infret;
                }
-       if (dd[0] + dd[1] == 0.) {
+       if (dval(&dd[0]) + dval(&dd[1]) == 0.) {
                b = buf;
 #ifndef IGNORE_ZERO_SIGN
                if (L[_0] & L[2+_0] & 0x80000000L)
@@ -99,16 +100,16 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
                return b;
                }
        if ((L[_0] & 0x7ff00000L) < (L[2+_0] & 0x7ff00000L)) {
-               ddx[1] = dd[0];
-               ddx[0] = dd[1];
+               dval(&ddx[1]) = dval(&dd[0]);
+               dval(&ddx[0]) = dval(&dd[1]);
                dd = ddx;
-               L = (ULong*)dd;
+               L = dd->L;
                }
-       z = d2b(dd[0], &ex, &bx);
-       if (dd[1] == 0.)
+       z = d2b(dval(&dd[0]), &ex, &bx);
+       if (dval(&dd[1]) == 0.)
                goto no_y;
        x = z;
-       y = d2b(dd[1], &ey, &by);
+       y = d2b(dval(&dd[1]), &ey, &by);
        if ( (i = ex - ey) !=0) {
                if (i > 0) {
                        x = lshift(x, i);
index 23d8b24..8367868 100644 (file)
@@ -88,6 +88,8 @@ g_dfmt(char *buf, double *d, int ndig, size_t bufsize)
        if (ndig <= 0)
                mode = 0;
        i = STRTOG_Normal;
+       if (sign)
+               i = STRTOG_Normal | STRTOG_Neg;
        s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
        return g__fmt(buf, s, se, decpt, sign, bufsize);
        }
index 3b64250..4686f63 100644 (file)
@@ -157,8 +157,9 @@ gdtoa
        int rdir, s2, s5, spec_case, try_quick;
        Long L;
        Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
-       double d, d2, ds, eps;
+       double d2, ds;
        char *s, *s0;
+       U d, eps;
 
 #ifndef MULTIPLE_THREADS
        if (dtoa_result) {
@@ -197,21 +198,21 @@ gdtoa
                return nrv_alloc("0", rve, 1);
                }
 
-       dval(d) = b2d(b, &i);
+       dval(&d) = b2d(b, &i);
        i = be + bbits - 1;
-       word0(d) &= Frac_mask1;
-       word0(d) |= Exp_11;
+       word0(&d) &= Frac_mask1;
+       word0(&d) |= Exp_11;
 #ifdef IBM
-       if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
-               dval(d) /= 1 << j;
+       if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
+               dval(&d) /= 1 << j;
 #endif
 
        /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
         * log10(x)      =  log(x) / log(10)
         *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
-        * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
+        * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
         *
-        * This suggests computing an approximation k to log10(d) by
+        * This suggests computing an approximation k to log10(&d) by
         *
         * k = (i - Bias)*0.301029995663981
         *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@@ -231,7 +232,7 @@ gdtoa
        i <<= 2;
        i += j;
 #endif
-       ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
+       ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 
        /* correct assumption about exponent range */
        if ((j = i) < 0)
@@ -246,13 +247,13 @@ gdtoa
 #ifdef IBM
        j = be + bbits - 1;
        if ( (j1 = j & 3) !=0)
-               dval(d) *= 1 << j1;
-       word0(d) += j << Exp_shift - 2 & Exp_mask;
+               dval(&d) *= 1 << j1;
+       word0(&d) += j << Exp_shift - 2 & Exp_mask;
 #else
-       word0(d) += (be + bbits - 1) << Exp_shift;
+       word0(&d) += (be + bbits - 1) << Exp_shift;
 #endif
        if (k >= 0 && k <= Ten_pmax) {
-               if (dval(d) < tens[k])
+               if (dval(&d) < tens[k])
                        k--;
                k_check = 0;
                }
@@ -283,10 +284,11 @@ gdtoa
                try_quick = 0;
                }
        leftright = 1;
+       ilim = ilim1 = -1;      /* Values for cases 0 and 1; done here to */
+                               /* silence erroneous "gcc -Wall" warning. */
        switch(mode) {
                case 0:
                case 1:
-                       ilim = ilim1 = -1;
                        i = (int)(nbits * .30103) + 3;
                        ndigits = 0;
                        break;
@@ -328,10 +330,10 @@ gdtoa
                /* Try to get by with floating-point arithmetic. */
 
                i = 0;
-               d2 = dval(d);
+               d2 = dval(&d);
 #ifdef IBM
-               if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
-                       dval(d) /= 1 << j;
+               if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
+                       dval(&d) /= 1 << j;
 #endif
                k0 = k;
                ilim0 = ilim;
@@ -342,7 +344,7 @@ gdtoa
                        if (j & Bletch) {
                                /* prevent overflows */
                                j &= Bletch - 1;
-                               dval(d) /= bigtens[n_bigtens-1];
+                               dval(&d) /= bigtens[n_bigtens-1];
                                ieps++;
                                }
                        for(; j; j >>= 1, i++)
@@ -354,30 +356,30 @@ gdtoa
                else  {
                        ds = 1.;
                        if ( (j1 = -k) !=0) {
-                               dval(d) *= tens[j1 & 0xf];
+                               dval(&d) *= tens[j1 & 0xf];
                                for(j = j1 >> 4; j; j >>= 1, i++)
                                        if (j & 1) {
                                                ieps++;
-                                               dval(d) *= bigtens[i];
+                                               dval(&d) *= bigtens[i];
                                                }
                                }
                        }
-               if (k_check && dval(d) < 1. && ilim > 0) {
+               if (k_check && dval(&d) < 1. && ilim > 0) {
                        if (ilim1 <= 0)
                                goto fast_failed;
                        ilim = ilim1;
                        k--;
-                       dval(d) *= 10.;
+                       dval(&d) *= 10.;
                        ieps++;
                        }
-               dval(eps) = ieps*dval(d) + 7.;
-               word0(eps) -= (P-1)*Exp_msk1;
+               dval(&eps) = ieps*dval(&d) + 7.;
+               word0(&eps) -= (P-1)*Exp_msk1;
                if (ilim == 0) {
                        S = mhi = 0;
-                       dval(d) -= 5.;
-                       if (dval(d) > dval(eps))
+                       dval(&d) -= 5.;
+                       if (dval(&d) > dval(&eps))
                                goto one_digit;
-                       if (dval(d) < -dval(eps))
+                       if (dval(&d) < -dval(&eps))
                                goto no_digits;
                        goto fast_failed;
                        }
@@ -386,38 +388,38 @@ gdtoa
                        /* Use Steele & White method of only
                         * generating digits needed.
                         */
-                       dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
+                       dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
                        for(i = 0;;) {
-                               L = (Long)(dval(d)/ds);
-                               dval(d) -= L*ds;
+                               L = (Long)(dval(&d)/ds);
+                               dval(&d) -= L*ds;
                                *s++ = '0' + (int)L;
-                               if (dval(d) < dval(eps)) {
-                                       if (dval(d))
+                               if (dval(&d) < dval(&eps)) {
+                                       if (dval(&d))
                                                inex = STRTOG_Inexlo;
                                        goto ret1;
                                        }
-                               if (ds - dval(d) < dval(eps))
+                               if (ds - dval(&d) < dval(&eps))
                                        goto bump_up;
                                if (++i >= ilim)
                                        break;
-                               dval(eps) *= 10.;
-                               dval(d) *= 10.;
+                               dval(&eps) *= 10.;
+                               dval(&d) *= 10.;
                                }
                        }
                else {
 #endif
                        /* Generate ilim digits, then fix them up. */
-                       dval(eps) *= tens[ilim-1];
-                       for(i = 1;; i++, dval(d) *= 10.) {
-                               if ( (L = (Long)(dval(d)/ds)) !=0)
-                                       dval(d) -= L*ds;
+                       dval(&eps) *= tens[ilim-1];
+                       for(i = 1;; i++, dval(&d) *= 10.) {
+                               if ( (L = (Long)(dval(&d)/ds)) !=0)
+                                       dval(&d) -= L*ds;
                                *s++ = '0' + (int)L;
                                if (i == ilim) {
                                        ds *= 0.5;
-                                       if (dval(d) > ds + dval(eps))
+                                       if (dval(&d) > ds + dval(&eps))
                                                goto bump_up;
-                                       else if (dval(d) < ds - dval(eps)) {
-                                               if (dval(d))
+                                       else if (dval(&d) < ds - dval(&eps)) {
+                                               if (dval(&d))
                                                        inex = STRTOG_Inexlo;
                                                goto clear_trailing0;
                                                }
@@ -429,7 +431,7 @@ gdtoa
 #endif
  fast_failed:
                s = s0;
-               dval(d) = d2;
+               dval(&d) = d2;
                k = k0;
                ilim = ilim0;
                }
@@ -441,22 +443,22 @@ gdtoa
                ds = tens[k];
                if (ndigits < 0 && ilim <= 0) {
                        S = mhi = 0;
-                       if (ilim < 0 || dval(d) <= 5*ds)
+                       if (ilim < 0 || dval(&d) <= 5*ds)
                                goto no_digits;
                        goto one_digit;
                        }
-               for(i = 1;; i++, dval(d) *= 10.) {
-                       L = dval(d) / ds;
-                       dval(d) -= L*ds;
+               for(i = 1;; i++, dval(&d) *= 10.) {
+                       L = dval(&d) / ds;
+                       dval(&d) -= L*ds;
 #ifdef Check_FLT_ROUNDS
                        /* If FLT_ROUNDS == 2, L will usually be high by 1 */
-                       if (dval(d) < 0) {
+                       if (dval(&d) < 0) {
                                L--;
-                               dval(d) += ds;
+                               dval(&d) += ds;
                                }
 #endif
                        *s++ = '0' + (int)L;
-                       if (dval(d) == 0.)
+                       if (dval(&d) == 0.)
                                break;
                        if (i == ilim) {
                                if (rdir) {
@@ -465,8 +467,13 @@ gdtoa
                                        inex = STRTOG_Inexlo;
                                        goto ret1;
                                        }
-                               dval(d) += dval(d);
-                               if (dval(d) > ds || dval(d) == ds && L & 1) {
+                               dval(&d) += dval(&d);
+#ifdef ROUND_BIASED
+                               if (dval(&d) >= ds)
+#else
+                               if (dval(&d) > ds || (dval(&d) == ds && L & 1))
+#endif
+                                       {
  bump_up:
                                        inex = STRTOG_Inexhi;
                                        while(*--s == '9')
@@ -560,28 +567,11 @@ gdtoa
         * and for all and pass them and a shift to quorem, so it
         * can do shifts and ors to compute the numerator for q.
         */
-#ifdef Pack_32
-       if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
-               i = 32 - i;
-#else
-       if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
-               i = 16 - i;
-#endif
-       if (i > 4) {
-               i -= 4;
-               b2 += i;
-               m2 += i;
-               s2 += i;
-               }
-       else if (i < 4) {
-               i += 28;
-               b2 += i;
-               m2 += i;
-               s2 += i;
-               }
-       if (b2 > 0)
+       i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
+       m2 += i;
+       if ((b2 += i) > 0)
                b = lshift(b, b2);
-       if (s2 > 0)
+       if ((s2 += i) > 0)
                S = lshift(S, s2);
        if (k_check) {
                if (cmp(b,S) < 0) {
@@ -646,11 +636,11 @@ gdtoa
                                goto ret;
                                }
 #endif
-                       if (j < 0 || j == 0 && !mode
+                       if (j < 0 || (j == 0 && !mode
 #ifndef ROUND_BIASED
                                                        && !(bits[0] & 1)
 #endif
-                                       ) {
+                                       )) {
                                if (rdir && (b->wds > 1 || b->x[0])) {
                                        if (rdir == 2) {
                                                inex = STRTOG_Inexlo;
@@ -673,7 +663,11 @@ gdtoa
                                if (j1 > 0) {
                                        b = lshift(b, 1);
                                        j1 = cmp(b, S);
-                                       if ((j1 > 0 || j1 == 0 && dig & 1)
+#ifdef ROUND_BIASED
+                                       if (j1 >= 0 /*)*/
+#else
+                                       if ((j1 > 0 || (j1 == 0 && dig & 1))
+#endif
                                        && dig++ == '9')
                                                goto round_9_up;
                                        inex = STRTOG_Inexhi;
@@ -718,13 +712,18 @@ gdtoa
        /* Round off last digit */
 
        if (rdir) {
-               if (rdir == 2 || b->wds <= 1 && !b->x[0])
+               if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
                        goto chopzeros;
                goto roundoff;
                }
        b = lshift(b, 1);
        j = cmp(b, S);
-       if (j > 0 || j == 0 && dig & 1) {
+#ifdef ROUND_BIASED
+       if (j >= 0)
+#else
+       if (j > 0 || (j == 0 && dig & 1))
+#endif
+               {
  roundoff:
                inex = STRTOG_Inexhi;
                while(*--s == '9')
index 3b038ff..2c188d6 100644 (file)
@@ -113,7 +113,12 @@ THIS SOFTWARE.
  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
  *     if memory is available and otherwise does something you deem
  *     appropriate.  If MALLOC is undefined, malloc will be invoked
- *     directly -- and assumed always to succeed.
+ *     directly -- and assumed always to succeed.  Similarly, if you
+ *     want something other than the system's free() to be called to
+ *     recycle memory acquired from MALLOC, #define FREE to be the
+ *     name of the alternate routine.  (FREE or free is only called in
+ *     pathological cases, e.g., in a gdtoa call after a gdtoa return in
+ *     mode 3 with thousands of digits requested.)
  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
  *     memory allocations from a private pool of memory when possible.
  *     When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
@@ -162,11 +167,6 @@ THIS SOFTWARE.
  * #define NO_STRING_H to use private versions of memcpy.
  *     On some K&R systems, it may also be necessary to
  *     #define DECLARE_SIZE_T in this case.
- * #define YES_ALIAS to permit aliasing certain double values with
- *     arrays of ULongs.  This leads to slightly better code with
- *     some compilers and was always used prior to 19990916, but it
- *     is not strictly legal and can cause trouble with aggressively
- *     optimizing compilers (e.g., gcc 2.95.1 under -O2).
  * #define USE_LOCALE to use the current locale's decimal_point value.
  */
 
@@ -270,25 +270,14 @@ Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
 
 typedef union { double d; ULong L[2]; } U;
 
-#ifdef YES_ALIAS
-#define dval(x) x
 #ifdef IEEE_8087
-#define word0(x) ((ULong *)&x)[1]
-#define word1(x) ((ULong *)&x)[0]
+#define word0(x) (x)->L[1]
+#define word1(x) (x)->L[0]
 #else
-#define word0(x) ((ULong *)&x)[0]
-#define word1(x) ((ULong *)&x)[1]
+#define word0(x) (x)->L[0]
+#define word1(x) (x)->L[1]
 #endif
-#else /* !YES_ALIAS */
-#ifdef IEEE_8087
-#define word0(x) ((U*)&x)->L[1]
-#define word1(x) ((U*)&x)->L[0]
-#else
-#define word0(x) ((U*)&x)->L[0]
-#define word1(x) ((U*)&x)->L[1]
-#endif
-#define dval(x) ((U*)&x)->d
-#endif /* YES_ALIAS */
+#define dval(x) (x)->d
 
 /* The following definition of Storeinc is appropriate for MIPS processors.
  * An alternative that might be better on some machines is
@@ -462,7 +451,7 @@ extern double rnd_prod(double, double), rnd_quot(double, double);
 #define FREE_DTOA_LOCK(n)      /*nothing*/
 #endif
 
-#define Kmax 15
+#define Kmax 9
 
  struct
 Bigint {
@@ -575,7 +564,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
  extern double strtod ANSI((const char *s00, char **se));
  extern Bigint *sum ANSI((Bigint*, Bigint*));
  extern int trailz ANSI((Bigint*));
- extern double ulp ANSI((double));
+ extern double ulp ANSI((U*));
 
 #ifdef __cplusplus
 }
index c8ac8ba..a9982c9 100644 (file)
@@ -57,7 +57,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
        static unsigned char *decimalpoint_cache;
        if (!(s0 = decimalpoint_cache)) {
                s0 = (unsigned char*)localeconv()->decimal_point;
-               if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
+               if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
                        strcpy(decimalpoint_cache, s0);
                        s0 = decimalpoint_cache;
                        }
@@ -199,7 +199,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
                return STRTOG_Normal | STRTOG_Inexlo;
                }
        n = s1 - s0 - 1;
-       for(k = 0; n > (1 << kshift-2) - 1; n >>= 1)
+       for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
                k++;
        b = Balloc(k);
        x = b->x;
@@ -314,7 +314,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
                        break;
                  case FPI_Round_near:
                        if (lostbits & 2
-                        && (lostbits & 1) | x[0] & 1)
+                        && (lostbits | x[0]) & 1)
                                up = 1;
                        break;
                  case FPI_Round_up:
@@ -333,8 +333,8 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
                                        irv =  STRTOG_Normal;
                                }
                        else if (b->wds > k
-                        || (n = nbits & kmask) !=0
-                            && hi0bits(x[k-1]) < 32-n) {
+                        || ((n = nbits & kmask) !=0
+                             && hi0bits(x[k-1]) < 32-n)) {
                                rshift(b,1);
                                if (++e > fpi->emax)
                                        goto ovfl;
index e07e069..a443721 100644 (file)
@@ -77,7 +77,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
        if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
         && *(CONST unsigned char*)(s+3) > ' ')
                s += 2;
-       while(c = *(CONST unsigned char*)++s) {
+       while((c = *(CONST unsigned char*)++s)) {
                if (!(h = hexdig[c])) {
                        if (c <= ' ') {
                                if (hd0 < havedig) {
@@ -109,7 +109,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
                                        *sp = s + 1;
                                        break;
                                        }
-                               } while(c = *++s);
+                               } while((c = *++s));
 #endif
                        return STRTOG_NaN;
                        }
@@ -120,7 +120,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
                        i = 1;
                        *--x = 0;
                        }
-               *x = (*x << 4) | h & 0xf;
+               *x = (*x << 4) | (h & 0xf);
                }
        if (!havedig)
                return STRTOG_NaN;
index 367eb9c..b1f18cd 100644 (file)
@@ -30,6 +30,8 @@ CFLAGS = -g
 .c.o:
        $(CC) -c $(CFLAGS) $*.c
 
+# invoke "make Printf" to add printf.o to gdtoa.a (if desired)
+
 all: arith.h gd_qnan.h gdtoa.a
 
 arith.h: arithchk.c
@@ -42,27 +44,36 @@ gd_qnan.h: arith.h qnan.c
        ./a.out >gd_qnan.h
        rm -f a.out qnan.o
 
-gdtoa.a: dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c g_ffmt.c\
-        g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c hexnan.c\
-        misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c\
-        strtoIx.c strtoIxL.c strtod.c strtodI.c strtodg.c strtof.c strtopQ.c\
-        strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c\
-        strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c sum.c ulp.c
+gdtoa.a:  dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
+        g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c\
+        hexnan.c misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c\
+        strtoIf.c strtoIg.c strtoIx.c strtoIxL.c strtod.c strtodI.c\
+        strtodg.c strtof.c strtopQ.c strtopd.c strtopdd.c strtopf.c\
+        strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c strtorf.c\
+        strtorx.c strtorxL.c sum.c ulp.c
        $(CC) -c $(CFLAGS) $?
        x=`echo $? | sed 's/\.c/.o/g'` && ar ruv gdtoa.a $$x && rm $$x
        ranlib gdtoa.a || true
 
+Printf: all printf.c
+       $(CC) -c $(CFLAGS) printf.c
+       ar ruv gdtoa.a printf.o
+       rm printf.o
+       touch Printf
+
 # If your system lacks ranlib, you do not need it.
 
-xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
-        g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h gdtoa_fltrnds.h gdtoaimp.h\
-        gethex.c gmisc.c hd_init.c hexnan.c makefile misc.c qnan.c smisc.c\
-        strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c strtoIxL.c\
-        strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c strtopQ.c strtopd.c\
-        strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c\
-        strtorf.c strtorx.c strtorxL.c sum.c ulp.c
+xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c\
+        g_dfmt.c g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h\
+        gdtoa_fltrnds.h gdtoaimp.h gethex.c gmisc.c hd_init.c hexnan.c\
+        makefile misc.c printf.c printf.c0 qnan.c smisc.c stdio1.h\
+        strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c\
+        strtoIxL.c strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c\
+        strtopQ.c strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c\
+        strtorQ.c strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c\
+        sum.c ulp.c
 
-# "make xsum.out" to check for transmission errors; source for xsum is
+# "make -r xsum.out" to check for transmission errors; source for xsum is
 # netlib's "xsum.c from f2c", e.g.,
 # ftp://netlib.bell-labs.com/netlib/f2c/xsum.c.gz
 
@@ -71,4 +82,4 @@ xsum.out: xsum0.out $(xs0)
        cmp xsum0.out xsum1.out && mv xsum1.out xsum.out || diff xsum[01].out
 
 clean:
-       rm -f arith.h gd_qnan.h *.[ao] xsum.out xsum1.out
+       rm -f arith.h gd_qnan.h *.[ao] Printf xsum.out xsum1.out
index b3ce7c9..e5f7b04 100644 (file)
@@ -55,7 +55,9 @@ Balloc
 #endif
 
        ACQUIRE_DTOA_LOCK(0);
-       if ( (rv = freelist[k]) !=0) {
+       /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
+       /* but this case seems very unlikely. */
+       if (k <= Kmax && (rv = freelist[k]) !=0) {
                freelist[k] = rv->next;
                }
        else {
@@ -65,7 +67,7 @@ Balloc
 #else
                len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
                        /sizeof(double);
-               if (pmem_next - private_mem + len <= PRIVATE_mem) {
+               if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
                        rv = (Bigint*)pmem_next;
                        pmem_next += len;
                        }
@@ -89,10 +91,18 @@ Bfree
 #endif
 {
        if (v) {
-               ACQUIRE_DTOA_LOCK(0);
-               v->next = freelist[v->k];
-               freelist[v->k] = v;
-               FREE_DTOA_LOCK(0);
+               if (v->k > Kmax)
+#ifdef FREE
+                       FREE((void*)v);
+#else
+                       free((void*)v);
+#endif
+               else {
+                       ACQUIRE_DTOA_LOCK(0);
+                       v->next = freelist[v->k];
+                       freelist[v->k] = v;
+                       FREE_DTOA_LOCK(0);
+                       }
                }
        }
 
@@ -104,8 +114,8 @@ lo0bits
        (ULong *y)
 #endif
 {
-       register int k;
-       register ULong x = *y;
+       int k;
+       ULong x = *y;
 
        if (x & 7) {
                if (x & 1)
@@ -204,12 +214,12 @@ multadd
  int
 hi0bits_D2A
 #ifdef KR_headers
-       (x) register ULong x;
+       (x) ULong x;
 #else
-       (register ULong x)
+       (ULong x)
 #endif
 {
-       register int k = 0;
+       int k = 0;
 
        if (!(x & 0xffff0000)) {
                k = 16;
@@ -612,12 +622,12 @@ b2d
 {
        ULong *xa, *xa0, w, y, z;
        int k;
-       double d;
+       U d;
 #ifdef VAX
        ULong d0, d1;
 #else
-#define d0 word0(d)
-#define d1 word1(d)
+#define d0 word0(&d)
+#define d1 word1(&d)
 #endif
 
        xa0 = a->x;
@@ -630,16 +640,16 @@ b2d
        *e = 32 - k;
 #ifdef Pack_32
        if (k < Ebits) {
-               d0 = Exp_1 | y >> Ebits - k;
+               d0 = Exp_1 | y >> (Ebits - k);
                w = xa > xa0 ? *--xa : 0;
-               d1 = y << (32-Ebits) + k | w >> Ebits - k;
+               d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
                goto ret_d;
                }
        z = xa > xa0 ? *--xa : 0;
        if (k -= Ebits) {
-               d0 = Exp_1 | y << k | z >> 32 - k;
+               d0 = Exp_1 | y << k | z >> (32 - k);
                y = xa > xa0 ? *--xa : 0;
-               d1 = z << k | y >> 32 - k;
+               d1 = z << k | y >> (32 - k);
                }
        else {
                d0 = Exp_1 | y;
@@ -663,10 +673,10 @@ b2d
 #endif
  ret_d:
 #ifdef VAX
-       word0(d) = d0 >> 16 | d0 << 16;
-       word1(d) = d1 >> 16 | d1 << 16;
+       word0(&d) = d0 >> 16 | d0 << 16;
+       word1(&d) = d1 >> 16 | d1 << 16;
 #endif
-       return dval(d);
+       return dval(&d);
        }
 #undef d0
 #undef d1
@@ -674,12 +684,13 @@ b2d
  Bigint *
 d2b
 #ifdef KR_headers
-       (d, e, bits) double d; int *e, *bits;
+       (dd, e, bits) double dd; int *e, *bits;
 #else
-       (double d, int *e, int *bits)
+       (double dd, int *e, int *bits)
 #endif
 {
        Bigint *b;
+       U d;
 #ifndef Sudden_Underflow
        int i;
 #endif
@@ -687,11 +698,14 @@ d2b
        ULong *x, y, z;
 #ifdef VAX
        ULong d0, d1;
-       d0 = word0(d) >> 16 | word0(d) << 16;
-       d1 = word1(d) >> 16 | word1(d) << 16;
 #else
-#define d0 word0(d)
-#define d1 word1(d)
+#define d0 word0(&d)
+#define d1 word1(&d)
+#endif
+       d.d = dd;
+#ifdef VAX
+       d0 = word0(&d) >> 16 | word0(&d) << 16;
+       d1 = word1(&d) >> 16 | word1(&d) << 16;
 #endif
 
 #ifdef Pack_32
@@ -715,7 +729,7 @@ d2b
 #ifdef Pack_32
        if ( (y = d1) !=0) {
                if ( (k = lo0bits(&y)) !=0) {
-                       x[0] = y | z << 32 - k;
+                       x[0] = y | z << (32 - k);
                        z >>= k;
                        }
                else
@@ -726,10 +740,6 @@ d2b
                     b->wds = (x[1] = z) !=0 ? 2 : 1;
                }
        else {
-#ifdef DEBUG
-               if (!z)
-                       Bug("Zero passed to d2b");
-#endif
                k = lo0bits(&z);
                x[0] = z;
 #ifndef Sudden_Underflow
@@ -788,7 +798,7 @@ d2b
 #endif
 #ifdef IBM
                *e = (de - Bias - (P-1) << 2) + k;
-               *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
+               *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
 #else
                *e = de - Bias - (P-1) + k;
                *bits = P - k;
@@ -841,7 +851,7 @@ strcp_D2A(a, b) char *a; char *b;
 strcp_D2A(char *a, CONST char *b)
 #endif
 {
-       while(*a = *b++)
+       while((*a = *b++))
                a++;
        return a;
        }
@@ -855,8 +865,8 @@ memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
 memcpy_D2A(void *a1, void *b1, size_t len)
 #endif
 {
-       register char *a = (char*)a1, *ae = a + len;
-       register char *b = (char*)b1, *a0 = a;
+       char *a = (char*)a1, *ae = a + len;
+       char *b = (char*)b1, *a0 = a;
        while(a < ae)
                *a++ = *b++;
        return a0;
index 50431fc..f4dbafb 100644 (file)
@@ -77,33 +77,33 @@ ratio
        (Bigint *a, Bigint *b)
 #endif
 {
-       double da, db;
+       U da, db;
        int k, ka, kb;
 
-       dval(da) = b2d(a, &ka);
-       dval(db) = b2d(b, &kb);
+       dval(&da) = b2d(a, &ka);
+       dval(&db) = b2d(b, &kb);
        k = ka - kb + ULbits*(a->wds - b->wds);
 #ifdef IBM
        if (k > 0) {
-               word0(da) += (k >> 2)*Exp_msk1;
+               word0(&da) += (k >> 2)*Exp_msk1;
                if (k &= 3)
-                       dval(da) *= 1 << k;
+                       dval(&da) *= 1 << k;
                }
        else {
                k = -k;
-               word0(db) += (k >> 2)*Exp_msk1;
+               word0(&db) += (k >> 2)*Exp_msk1;
                if (k &= 3)
-                       dval(db) *= 1 << k;
+                       dval(&db) *= 1 << k;
                }
 #else
        if (k > 0)
-               word0(da) += k*Exp_msk1;
+               word0(&da) += k*Exp_msk1;
        else {
                k = -k;
-               word0(db) += k*Exp_msk1;
+               word0(&db) += k*Exp_msk1;
                }
 #endif
-       return dval(da) / dval(db);
+       return dval(&da) / dval(&db);
        }
 
 #ifdef INFNAN_CHECK
index 90c88db..6a17760 100644 (file)
@@ -74,7 +74,7 @@ strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
                        goto swapcheck;
                        }
                if (b1->wds > nw
-                || nb1 && b1->x[nw1] & 1L << nb1) {
+                || (nb1 && b1->x[nw1] & 1L << nb1)) {
                        if (++e1 > fpi->emax)
                                rv1 = STRTOG_Infinite | STRTOG_Inexhi;
                        rshift(b1, 1);
index b8287eb..9783371 100644 (file)
@@ -57,6 +57,28 @@ static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
 #define Rounding Flt_Rounds
 #endif
 
+#ifdef Avoid_Underflow /*{*/
+ static double
+sulp
+#ifdef KR_headers
+       (x, scale) U *x; int scale;
+#else
+       (U *x, int scale)
+#endif
+{
+       U u;
+       double rv;
+       int i;
+
+       rv = ulp(x);
+       if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
+               return rv; /* Is there an example where i <= 0 ? */
+       word0(&u) = Exp_1 + (i << Exp_shift);
+       word1(&u) = 0;
+       return rv * u.d;
+       }
+#endif /*}*/
+
  double
 strtod
 #ifdef KR_headers
@@ -71,10 +93,14 @@ strtod
        int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
                 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
        CONST char *s, *s0, *s1;
-       double aadj, aadj1, adj, rv, rv0;
+       double aadj;
        Long L;
+       U adj, aadj1, rv, rv0;
        ULong y, z;
        Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
+#ifdef Avoid_Underflow
+       ULong Lsb, Lsb1;
+#endif
 #ifdef SET_INEXACT
        int inexact, oldinexact;
 #endif
@@ -88,7 +114,7 @@ strtod
        static int dplen;
        if (!(s0 = decimalpoint_cache)) {
                s0 = localeconv()->decimal_point;
-               if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
+               if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
                        strcpy(decimalpoint_cache, s0);
                        s0 = decimalpoint_cache;
                        }
@@ -115,7 +141,7 @@ strtod
 #endif /*}*/
 
        sign = nz0 = nz = decpt = 0;
-       dval(rv) = 0.;
+       dval(&rv) = 0.;
        for(s = s00;;s++) switch(*s) {
                case '-':
                        sign = 1;
@@ -147,20 +173,12 @@ strtod
                  case 'x':
                  case 'X':
                        {
-#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/
+#ifdef Honor_FLT_ROUNDS
                        FPI fpi1 = fpi;
-#ifdef Honor_FLT_ROUNDS /*{{*/
                        fpi1.rounding = Rounding;
-#else /*}{*/
-                       switch(fegetround()) {
-                         case FE_TOWARDZERO:   fpi1.rounding = 0; break;
-                         case FE_UPWARD:       fpi1.rounding = 2; break;
-                         case FE_DOWNWARD:     fpi1.rounding = 3;
-                         }
-#endif /*}}*/
-#else /*}{*/
+#else
 #define fpi1 fpi
-#endif /*}}*/
+#endif
                        switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
                          case STRTOG_NoNumber:
                                s = s00;
@@ -285,8 +303,8 @@ strtod
                                        --s;
                                        if (!match(&s,"inity"))
                                                ++s;
-                                       word0(rv) = 0x7ff00000;
-                                       word1(rv) = 0;
+                                       word0(&rv) = 0x7ff00000;
+                                       word1(&rv) = 0;
                                        goto ret;
                                        }
                                break;
@@ -297,13 +315,13 @@ strtod
                                        if (*s == '(' /*)*/
                                         && hexnan(&s, &fpinan, bits)
                                                        == STRTOG_NaNbits) {
-                                               word0(rv) = 0x7ff00000 | bits[1];
-                                               word1(rv) = bits[0];
+                                               word0(&rv) = 0x7ff00000 | bits[1];
+                                               word1(&rv) = bits[0];
                                                }
                                        else {
 #endif
-                                               word0(rv) = NAN_WORD0;
-                                               word1(rv) = NAN_WORD1;
+                                               word0(&rv) = NAN_WORD0;
+                                               word1(&rv) = NAN_WORD1;
 #ifndef No_Hex_NaN
                                                }
 #endif
@@ -327,13 +345,13 @@ strtod
        if (!nd0)
                nd0 = nd;
        k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-       dval(rv) = y;
+       dval(&rv) = y;
        if (k > 9) {
 #ifdef SET_INEXACT
                if (k > DBL_DIG)
                        oldinexact = get_inexact();
 #endif
-               dval(rv) = tens[k - 9] * dval(rv) + z;
+               dval(&rv) = tens[k - 9] * dval(&rv) + z;
                }
        bd0 = 0;
        if (nd <= DBL_DIG
@@ -353,11 +371,11 @@ strtod
 #ifdef Honor_FLT_ROUNDS
                                /* round correctly FLT_ROUNDS = 2 or 3 */
                                if (sign) {
-                                       rv = -rv;
+                                       rv.d = -rv.d;
                                        sign = 0;
                                        }
 #endif
-                               /* rv = */ rounded_product(dval(rv), tens[e]);
+                               /* rv = */ rounded_product(dval(&rv), tens[e]);
                                goto ret;
 #endif
                                }
@@ -369,25 +387,25 @@ strtod
 #ifdef Honor_FLT_ROUNDS
                                /* round correctly FLT_ROUNDS = 2 or 3 */
                                if (sign) {
-                                       rv = -rv;
+                                       rv.d = -rv.d;
                                        sign = 0;
                                        }
 #endif
                                e -= i;
-                               dval(rv) *= tens[i];
+                               dval(&rv) *= tens[i];
 #ifdef VAX
                                /* VAX exponent range is so narrow we must
                                 * worry about overflow here...
                                 */
  vax_ovfl_check:
-                               word0(rv) -= P*Exp_msk1;
-                               /* rv = */ rounded_product(dval(rv), tens[e]);
-                               if ((word0(rv) & Exp_mask)
+                               word0(&rv) -= P*Exp_msk1;
+                               /* rv = */ rounded_product(dval(&rv), tens[e]);
+                               if ((word0(&rv) & Exp_mask)
                                 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
                                        goto ovfl;
-                               word0(rv) += P*Exp_msk1;
+                               word0(&rv) += P*Exp_msk1;
 #else
-                               /* rv = */ rounded_product(dval(rv), tens[e]);
+                               /* rv = */ rounded_product(dval(&rv), tens[e]);
 #endif
                                goto ret;
                                }
@@ -397,11 +415,11 @@ strtod
 #ifdef Honor_FLT_ROUNDS
                        /* round correctly FLT_ROUNDS = 2 or 3 */
                        if (sign) {
-                               rv = -rv;
+                               rv.d = -rv.d;
                                sign = 0;
                                }
 #endif
-                       /* rv = */ rounded_quotient(dval(rv), tens[-e]);
+                       /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
                        goto ret;
                        }
 #endif
@@ -432,67 +450,73 @@ strtod
 
        if (e1 > 0) {
                if ( (i = e1 & 15) !=0)
-                       dval(rv) *= tens[i];
+                       dval(&rv) *= tens[i];
                if (e1 &= ~15) {
                        if (e1 > DBL_MAX_10_EXP) {
  ovfl:
-#ifndef NO_ERRNO
-                               errno = ERANGE;
-#endif
                                /* Can't trust HUGE_VAL */
 #ifdef IEEE_Arith
 #ifdef Honor_FLT_ROUNDS
                                switch(Rounding) {
                                  case 0: /* toward 0 */
                                  case 3: /* toward -infinity */
-                                       word0(rv) = Big0;
-                                       word1(rv) = Big1;
+                                       word0(&rv) = Big0;
+                                       word1(&rv) = Big1;
                                        break;
                                  default:
-                                       word0(rv) = Exp_mask;
-                                       word1(rv) = 0;
+                                       word0(&rv) = Exp_mask;
+                                       word1(&rv) = 0;
                                  }
 #else /*Honor_FLT_ROUNDS*/
-                               word0(rv) = Exp_mask;
-                               word1(rv) = 0;
+                               word0(&rv) = Exp_mask;
+                               word1(&rv) = 0;
 #endif /*Honor_FLT_ROUNDS*/
 #ifdef SET_INEXACT
                                /* set overflow bit */
-                               dval(rv0) = 1e300;
-                               dval(rv0) *= dval(rv0);
+                               dval(&rv0) = 1e300;
+                               dval(&rv0) *= dval(&rv0);
 #endif
 #else /*IEEE_Arith*/
-                               word0(rv) = Big0;
-                               word1(rv) = Big1;
+                               word0(&rv) = Big0;
+                               word1(&rv) = Big1;
 #endif /*IEEE_Arith*/
-                               if (bd0)
-                                       goto retfree;
+ range_err:
+                               if (bd0) {
+                                       Bfree(bb);
+                                       Bfree(bd);
+                                       Bfree(bs);
+                                       Bfree(bd0);
+                                       Bfree(delta);
+                                       }
+#ifndef NO_ERRNO
+                               errno = ERANGE;
+#endif
                                goto ret;
                                }
                        e1 >>= 4;
                        for(j = 0; e1 > 1; j++, e1 >>= 1)
                                if (e1 & 1)
-                                       dval(rv) *= bigtens[j];
+                                       dval(&rv) *= bigtens[j];
                /* The last multiplication could overflow. */
-                       word0(rv) -= P*Exp_msk1;
-                       dval(rv) *= bigtens[j];
-                       if ((z = word0(rv) & Exp_mask)
+                       word0(&rv) -= P*Exp_msk1;
+                       dval(&rv) *= bigtens[j];
+                       if ((z = word0(&rv) & Exp_mask)
                         > Exp_msk1*(DBL_MAX_EXP+Bias-P))
                                goto ovfl;
                        if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
                                /* set to largest number */
                                /* (Can't trust DBL_MAX) */
-                               word0(rv) = Big0;
-                               word1(rv) = Big1;
+                               word0(&rv) = Big0;
+                               word1(&rv) = Big1;
                                }
                        else
-                               word0(rv) += P*Exp_msk1;
+                               word0(&rv) += P*Exp_msk1;
                        }
                }
        else if (e1 < 0) {
                e1 = -e1;
                if ( (i = e1 & 15) !=0)
-                       dval(rv) /= tens[i];
+                       dval(&rv) /= tens[i];
                if (e1 >>= 4) {
                        if (e1 >= 1 << n_bigtens)
                                goto undfl;
@@ -501,44 +525,39 @@ strtod
                                scale = 2*P;
                        for(j = 0; e1 > 0; j++, e1 >>= 1)
                                if (e1 & 1)
-                                       dval(rv) *= tinytens[j];
-                       if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
+                                       dval(&rv) *= tinytens[j];
+                       if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
                                                >> Exp_shift)) > 0) {
                                /* scaled rv is denormal; zap j low bits */
                                if (j >= 32) {
-                                       word1(rv) = 0;
+                                       word1(&rv) = 0;
                                        if (j >= 53)
-                                        word0(rv) = (P+2)*Exp_msk1;
+                                        word0(&rv) = (P+2)*Exp_msk1;
                                        else
-                                        word0(rv) &= 0xffffffff << j-32;
+                                        word0(&rv) &= 0xffffffff << (j-32);
                                        }
                                else
-                                       word1(rv) &= 0xffffffff << j;
+                                       word1(&rv) &= 0xffffffff << j;
                                }
 #else
                        for(j = 0; e1 > 1; j++, e1 >>= 1)
                                if (e1 & 1)
-                                       dval(rv) *= tinytens[j];
+                                       dval(&rv) *= tinytens[j];
                        /* The last multiplication could underflow. */
-                       dval(rv0) = dval(rv);
-                       dval(rv) *= tinytens[j];
-                       if (!dval(rv)) {
-                               dval(rv) = 2.*dval(rv0);
-                               dval(rv) *= tinytens[j];
+                       dval(&rv0) = dval(&rv);
+                       dval(&rv) *= tinytens[j];
+                       if (!dval(&rv)) {
+                               dval(&rv) = 2.*dval(&rv0);
+                               dval(&rv) *= tinytens[j];
 #endif
-                               if (!dval(rv)) {
+                               if (!dval(&rv)) {
  undfl:
-                                       dval(rv) = 0.;
-#ifndef NO_ERRNO
-                                       errno = ERANGE;
-#endif
-                                       if (bd0)
-                                               goto retfree;
-                                       goto ret;
+                                       dval(&rv) = 0.;
+                                       goto range_err;
                                        }
 #ifndef Avoid_Underflow
-                               word0(rv) = Tiny0;
-                               word1(rv) = Tiny1;
+                               word0(&rv) = Tiny0;
+                               word1(&rv) = Tiny1;
                                /* The refinement below will clean
                                 * this approximation up.
                                 */
@@ -556,7 +575,7 @@ strtod
        for(;;) {
                bd = Balloc(bd0->k);
                Bcopy(bd, bd0);
-               bb = d2b(dval(rv), &bbe, &bbbits);      /* rv = bb * 2^bbe */
+               bb = d2b(dval(&rv), &bbe, &bbbits);     /* rv = bb * 2^bbe */
                bs = i2b(1);
 
                if (e >= 0) {
@@ -577,12 +596,19 @@ strtod
                        bs2++;
 #endif
 #ifdef Avoid_Underflow
+               Lsb = LSB;
+               Lsb1 = 0;
                j = bbe - scale;
                i = j + bbbits - 1;     /* logb(rv) */
-               if (i < Emin)   /* denormal */
-                       j += P - Emin;
-               else
-                       j = P + 1 - bbbits;
+               j = P + 1 - bbbits;
+               if (i < Emin) { /* denormal */
+                       i = Emin - i;
+                       j -= i;
+                       if (i < 32)
+                               Lsb <<= i;
+                       else
+                               Lsb1 = Lsb << (i-32);
+                       }
 #else /*Avoid_Underflow*/
 #ifdef Sudden_Underflow
 #ifdef IBM
@@ -592,7 +618,7 @@ strtod
 #endif
 #else /*Sudden_Underflow*/
                j = bbe;
-               i = j + bbbits - 1;     /* logb(rv) */
+               i = j + bbbits - 1;     /* logb(&rv) */
                if (i < Emin)   /* denormal */
                        j += P - Emin;
                else
@@ -643,15 +669,15 @@ strtod
                                        }
                                if (Rounding) {
                                        if (dsign) {
-                                               adj = 1.;
+                                               dval(&adj) = 1.;
                                                goto apply_adj;
                                                }
                                        }
                                else if (!dsign) {
-                                       adj = -1.;
-                                       if (!word1(rv)
-                                        && !(word0(rv) & Frac_mask)) {
-                                               y = word0(rv) & Exp_mask;
+                                       dval(&adj) = -1.;
+                                       if (!word1(&rv)
+                                        && !(word0(&rv) & Frac_mask)) {
+                                               y = word0(&rv) & Exp_mask;
 #ifdef Avoid_Underflow
                                                if (!scale || y > 2*P*Exp_msk1)
 #else
@@ -660,66 +686,66 @@ strtod
                                                  {
                                                  delta = lshift(delta,Log2P);
                                                  if (cmp(delta, bs) <= 0)
-                                                       adj = -0.5;
+                                                       dval(&adj) = -0.5;
                                                  }
                                                }
  apply_adj:
 #ifdef Avoid_Underflow
-                                       if (scale && (y = word0(rv) & Exp_mask)
+                                       if (scale && (y = word0(&rv) & Exp_mask)
                                                <= 2*P*Exp_msk1)
-                                         word0(adj) += (2*P+1)*Exp_msk1 - y;
+                                         word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
-                                       if ((word0(rv) & Exp_mask) <=
+                                       if ((word0(&rv) & Exp_mask) <=
                                                        P*Exp_msk1) {
-                                               word0(rv) += P*Exp_msk1;
-                                               dval(rv) += adj*ulp(dval(rv));
-                                               word0(rv) -= P*Exp_msk1;
+                                               word0(&rv) += P*Exp_msk1;
+                                               dval(&rv) += adj*ulp(&rv);
+                                               word0(&rv) -= P*Exp_msk1;
                                                }
                                        else
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
-                                       dval(rv) += adj*ulp(dval(rv));
+                                       dval(&rv) += adj.d*ulp(&rv);
                                        }
                                break;
                                }
-                       adj = ratio(delta, bs);
-                       if (adj < 1.)
-                               adj = 1.;
-                       if (adj <= 0x7ffffffe) {
-                               /* adj = Rounding ? ceil(adj) : floor(adj); */
-                               y = adj;
-                               if (y != adj) {
+                       dval(&adj) = ratio(delta, bs);
+                       if (adj.d < 1.)
+                               dval(&adj) = 1.;
+                       if (adj.d <= 0x7ffffffe) {
+                               /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
+                               y = adj.d;
+                               if (y != adj.d) {
                                        if (!((Rounding>>1) ^ dsign))
                                                y++;
-                                       adj = y;
+                                       dval(&adj) = y;
                                        }
                                }
 #ifdef Avoid_Underflow
-                       if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
-                               word0(adj) += (2*P+1)*Exp_msk1 - y;
+                       if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
+                               word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
-                       if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
-                               word0(rv) += P*Exp_msk1;
-                               adj *= ulp(dval(rv));
+                       if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
+                               word0(&rv) += P*Exp_msk1;
+                               dval(&adj) *= ulp(&rv);
                                if (dsign)
-                                       dval(rv) += adj;
+                                       dval(&rv) += adj;
                                else
-                                       dval(rv) -= adj;
-                               word0(rv) -= P*Exp_msk1;
+                                       dval(&rv) -= adj;
+                               word0(&rv) -= P*Exp_msk1;
                                goto cont;
                                }
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
-                       adj *= ulp(dval(rv));
+                       dval(&adj) *= ulp(&rv);
                        if (dsign) {
-                               if (word0(rv) == Big0 && word1(rv) == Big1)
+                               if (word0(&rv) == Big0 && word1(&rv) == Big1)
                                        goto ovfl;
-                               dval(rv) += adj;
+                               dval(&rv) += adj.d;
                                }
                        else
-                               dval(rv) -= adj;
+                               dval(&rv) -= adj.d;
                        goto cont;
                        }
 #endif /*Honor_FLT_ROUNDS*/
@@ -728,12 +754,12 @@ strtod
                        /* Error is less than half an ulp -- check for
                         * special case of mantissa a power of two.
                         */
-                       if (dsign || word1(rv) || word0(rv) & Bndry_mask
+                       if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
 #ifdef IEEE_Arith
 #ifdef Avoid_Underflow
-                        || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
+                        || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
 #else
-                        || (word0(rv) & Exp_mask) <= Exp_msk1
+                        || (word0(&rv) & Exp_mask) <= Exp_msk1
 #endif
 #endif
                                ) {
@@ -758,32 +784,34 @@ strtod
                if (i == 0) {
                        /* exactly half-way between */
                        if (dsign) {
-                               if ((word0(rv) & Bndry_mask1) == Bndry_mask1
-                                &&  word1(rv) == (
+                               if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
+                                &&  word1(&rv) == (
 #ifdef Avoid_Underflow
-                       (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
+                       (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
                ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
 #endif
                                                   0xffffffff)) {
                                        /*boundary case -- increment exponent*/
-                                       word0(rv) = (word0(rv) & Exp_mask)
+                                       if (word0(&rv) == Big0 && word1(&rv) == Big1)
+                                               goto ovfl;
+                                       word0(&rv) = (word0(&rv) & Exp_mask)
                                                + Exp_msk1
 #ifdef IBM
                                                | Exp_msk1 >> 4
 #endif
                                                ;
-                                       word1(rv) = 0;
+                                       word1(&rv) = 0;
 #ifdef Avoid_Underflow
                                        dsign = 0;
 #endif
                                        break;
                                        }
                                }
-                       else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
+                       else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
  drop_down:
                                /* boundary case -- decrement exponent */
 #ifdef Sudden_Underflow /*{{*/
-                               L = word0(rv) & Exp_mask;
+                               L = word0(&rv) & Exp_mask;
 #ifdef IBM
                                if (L <  Exp_msk1)
 #else
@@ -798,7 +826,7 @@ strtod
 #else /*Sudden_Underflow}{*/
 #ifdef Avoid_Underflow
                                if (scale) {
-                                       L = word0(rv) & Exp_mask;
+                                       L = word0(&rv) & Exp_mask;
                                        if (L <= (2*P+1)*Exp_msk1) {
                                                if (L > (P+2)*Exp_msk1)
                                                        /* round even ==> */
@@ -809,10 +837,10 @@ strtod
                                                }
                                        }
 #endif /*Avoid_Underflow*/
-                               L = (word0(rv) & Exp_mask) - Exp_msk1;
+                               L = (word0(&rv) & Exp_mask) - Exp_msk1;
 #endif /*Sudden_Underflow}}*/
-                               word0(rv) = L | Bndry_mask1;
-                               word1(rv) = 0xffffffff;
+                               word0(&rv) = L | Bndry_mask1;
+                               word1(&rv) = 0xffffffff;
 #ifdef IBM
                                goto cont;
 #else
@@ -820,16 +848,33 @@ strtod
 #endif
                                }
 #ifndef ROUND_BIASED
-                       if (!(word1(rv) & LSB))
+#ifdef Avoid_Underflow
+                       if (Lsb1) {
+                               if (!(word0(&rv) & Lsb1))
+                                       break;
+                               }
+                       else if (!(word1(&rv) & Lsb))
+                               break;
+#else
+                       if (!(word1(&rv) & LSB))
                                break;
+#endif
 #endif
                        if (dsign)
-                               dval(rv) += ulp(dval(rv));
+#ifdef Avoid_Underflow
+                               dval(&rv) += sulp(&rv, scale);
+#else
+                               dval(&rv) += ulp(&rv);
+#endif
 #ifndef ROUND_BIASED
                        else {
-                               dval(rv) -= ulp(dval(rv));
+#ifdef Avoid_Underflow
+                               dval(&rv) -= sulp(&rv, scale);
+#else
+                               dval(&rv) -= ulp(&rv);
+#endif
 #ifndef Sudden_Underflow
-                               if (!dval(rv))
+                               if (!dval(&rv))
                                        goto undfl;
 #endif
                                }
@@ -841,14 +886,14 @@ strtod
                        }
                if ((aadj = ratio(delta, bs)) <= 2.) {
                        if (dsign)
-                               aadj = aadj1 = 1.;
-                       else if (word1(rv) || word0(rv) & Bndry_mask) {
+                               aadj = dval(&aadj1) = 1.;
+                       else if (word1(&rv) || word0(&rv) & Bndry_mask) {
 #ifndef Sudden_Underflow
-                               if (word1(rv) == Tiny1 && !word0(rv))
+                               if (word1(&rv) == Tiny1 && !word0(&rv))
                                        goto undfl;
 #endif
                                aadj = 1.;
-                               aadj1 = -1.;
+                               dval(&aadj1) = -1.;
                                }
                        else {
                                /* special case -- power of FLT_RADIX to be */
@@ -858,45 +903,45 @@ strtod
                                        aadj = 1./FLT_RADIX;
                                else
                                        aadj *= 0.5;
-                               aadj1 = -aadj;
+                               dval(&aadj1) = -aadj;
                                }
                        }
                else {
                        aadj *= 0.5;
-                       aadj1 = dsign ? aadj : -aadj;
+                       dval(&aadj1) = dsign ? aadj : -aadj;
 #ifdef Check_FLT_ROUNDS
                        switch(Rounding) {
                                case 2: /* towards +infinity */
-                                       aadj1 -= 0.5;
+                                       dval(&aadj1) -= 0.5;
                                        break;
                                case 0: /* towards 0 */
                                case 3: /* towards -infinity */
-                                       aadj1 += 0.5;
+                                       dval(&aadj1) += 0.5;
                                }
 #else
                        if (Flt_Rounds == 0)
-                               aadj1 += 0.5;
+                               dval(&aadj1) += 0.5;
 #endif /*Check_FLT_ROUNDS*/
                        }
-               y = word0(rv) & Exp_mask;
+               y = word0(&rv) & Exp_mask;
 
                /* Check for overflow */
 
                if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
-                       dval(rv0) = dval(rv);
-                       word0(rv) -= P*Exp_msk1;
-                       adj = aadj1 * ulp(dval(rv));
-                       dval(rv) += adj;
-                       if ((word0(rv) & Exp_mask) >=
+                       dval(&rv0) = dval(&rv);
+                       word0(&rv) -= P*Exp_msk1;
+                       dval(&adj) = dval(&aadj1) * ulp(&rv);
+                       dval(&rv) += dval(&adj);
+                       if ((word0(&rv) & Exp_mask) >=
                                        Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
-                               if (word0(rv0) == Big0 && word1(rv0) == Big1)
+                               if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
                                        goto ovfl;
-                               word0(rv) = Big0;
-                               word1(rv) = Big1;
+                               word0(&rv) = Big0;
+                               word1(&rv) = Big1;
                                goto cont;
                                }
                        else
-                               word0(rv) += P*Exp_msk1;
+                               word0(&rv) += P*Exp_msk1;
                        }
                else {
 #ifdef Avoid_Underflow
@@ -905,58 +950,58 @@ strtod
                                        if ((z = aadj) <= 0)
                                                z = 1;
                                        aadj = z;
-                                       aadj1 = dsign ? aadj : -aadj;
+                                       dval(&aadj1) = dsign ? aadj : -aadj;
                                        }
-                               word0(aadj1) += (2*P+1)*Exp_msk1 - y;
+                               word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
                                }
-                       adj = aadj1 * ulp(dval(rv));
-                       dval(rv) += adj;
+                       dval(&adj) = dval(&aadj1) * ulp(&rv);
+                       dval(&rv) += dval(&adj);
 #else
 #ifdef Sudden_Underflow
-                       if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
-                               dval(rv0) = dval(rv);
-                               word0(rv) += P*Exp_msk1;
-                               adj = aadj1 * ulp(dval(rv));
-                               dval(rv) += adj;
+                       if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
+                               dval(&rv0) = dval(&rv);
+                               word0(&rv) += P*Exp_msk1;
+                               dval(&adj) = dval(&aadj1) * ulp(&rv);
+                               dval(&rv) += adj;
 #ifdef IBM
-                               if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
+                               if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
 #else
-                               if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
+                               if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
 #endif
                                        {
-                                       if (word0(rv0) == Tiny0
-                                        && word1(rv0) == Tiny1)
+                                       if (word0(&rv0) == Tiny0
+                                        && word1(&rv0) == Tiny1)
                                                goto undfl;
-                                       word0(rv) = Tiny0;
-                                       word1(rv) = Tiny1;
+                                       word0(&rv) = Tiny0;
+                                       word1(&rv) = Tiny1;
                                        goto cont;
                                        }
                                else
-                                       word0(rv) -= P*Exp_msk1;
+                                       word0(&rv) -= P*Exp_msk1;
                                }
                        else {
-                               adj = aadj1 * ulp(dval(rv));
-                               dval(rv) += adj;
+                               dval(&adj) = dval(&aadj1) * ulp(&rv);
+                               dval(&rv) += adj;
                                }
 #else /*Sudden_Underflow*/
-                       /* Compute adj so that the IEEE rounding rules will
-                        * correctly round rv + adj in some half-way cases.
-                        * If rv * ulp(rv) is denormalized (i.e.,
+                       /* Compute dval(&adj) so that the IEEE rounding rules will
+                        * correctly round rv + dval(&adj) in some half-way cases.
+                        * If rv * ulp(&rv) is denormalized (i.e.,
                         * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
                         * trouble from bits lost to denormalization;
                         * example: 1.2e-307 .
                         */
                        if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
-                               aadj1 = (double)(int)(aadj + 0.5);
+                               dval(&aadj1) = (double)(int)(aadj + 0.5);
                                if (!dsign)
-                                       aadj1 = -aadj1;
+                                       dval(&aadj1) = -dval(&aadj1);
                                }
-                       adj = aadj1 * ulp(dval(rv));
-                       dval(rv) += adj;
+                       dval(&adj) = dval(&aadj1) * ulp(&rv);
+                       dval(&rv) += adj;
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
                        }
-               z = word0(rv) & Exp_mask;
+               z = word0(&rv) & Exp_mask;
 #ifndef SET_INEXACT
 #ifdef Avoid_Underflow
                if (!scale)
@@ -966,7 +1011,7 @@ strtod
                        L = (Long)aadj;
                        aadj -= L;
                        /* The tolerances below are conservative. */
-                       if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
+                       if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
                                if (aadj < .4999999 || aadj > .5000001)
                                        break;
                                }
@@ -980,12 +1025,17 @@ strtod
                Bfree(bs);
                Bfree(delta);
                }
+       Bfree(bb);
+       Bfree(bd);
+       Bfree(bs);
+       Bfree(bd0);
+       Bfree(delta);
 #ifdef SET_INEXACT
        if (inexact) {
                if (!oldinexact) {
-                       word0(rv0) = Exp_1 + (70 << Exp_shift);
-                       word1(rv0) = 0;
-                       dval(rv0) += 1.;
+                       word0(&rv0) = Exp_1 + (70 << Exp_shift);
+                       word1(&rv0) = 0;
+                       dval(&rv0) += 1.;
                        }
                }
        else if (!oldinexact)
@@ -993,36 +1043,30 @@ strtod
 #endif
 #ifdef Avoid_Underflow
        if (scale) {
-               word0(rv0) = Exp_1 - 2*P*Exp_msk1;
-               word1(rv0) = 0;
-               dval(rv) *= dval(rv0);
+               word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
+               word1(&rv0) = 0;
+               dval(&rv) *= dval(&rv0);
 #ifndef NO_ERRNO
                /* try to avoid the bug of testing an 8087 register value */
 #ifdef IEEE_Arith
-               if (!(word0(rv) & Exp_mask))
+               if (!(word0(&rv) & Exp_mask))
 #else
-               if (word0(rv) == 0 && word1(rv) == 0)
+               if (word0(&rv) == 0 && word1(&rv) == 0)
 #endif
                        errno = ERANGE;
 #endif
                }
 #endif /* Avoid_Underflow */
 #ifdef SET_INEXACT
-       if (inexact && !(word0(rv) & Exp_mask)) {
+       if (inexact && !(word0(&rv) & Exp_mask)) {
                /* set underflow bit */
-               dval(rv0) = 1e-300;
-               dval(rv0) *= dval(rv0);
+               dval(&rv0) = 1e-300;
+               dval(&rv0) *= dval(&rv0);
                }
 #endif
- retfree:
-       Bfree(bb);
-       Bfree(bd);
-       Bfree(bs);
-       Bfree(bd0);
-       Bfree(delta);
  ret:
        if (se)
                *se = (char *)s;
-       return sign ? -dval(rv) : dval(rv);
+       return sign ? -dval(&rv) : dval(&rv);
        }
 
index 98f8891..0b7b8a4 100644 (file)
@@ -33,16 +33,16 @@ THIS SOFTWARE.
 
  static double
 #ifdef KR_headers
-ulpdown(d) double *d;
+ulpdown(d) U *d;
 #else
-ulpdown(double *d)
+ulpdown(U *d)
 #endif
 {
        double u;
-       ULong *L = (ULong*)d;
+       ULong *L = d->L;
 
-       u = ulp(*d);
-       if (!(L[_1] | L[_0] & 0xfffff)
+       u = ulp(d);
+       if (!(L[_1] | (L[_0] & 0xfffff))
         && (L[_0] & 0x7ff00000) > 0x00100000)
                u *= 0.5;
        return u;
@@ -59,10 +59,6 @@ strtodI(CONST char *s, char **sp, double *dd)
        ULong bits[2], sign;
        Long exp;
        int j, k;
-       typedef union {
-               double d[2];
-               ULong L[4];
-               } U;
        U *u;
 
        k = strtodg(s, sp, &fpi, &exp, bits);
@@ -70,17 +66,17 @@ strtodI(CONST char *s, char **sp, double *dd)
        sign = k & STRTOG_Neg ? 0x80000000L : 0;
        switch(k & STRTOG_Retmask) {
          case STRTOG_NoNumber:
-               u->d[0] = u->d[1] = 0.;
+               dval(&u[0]) = dval(&u[1]) = 0.;
                break;
 
          case STRTOG_Zero:
-               u->d[0] = u->d[1] = 0.;
+               dval(&u[0]) = dval(&u[1]) = 0.;
 #ifdef Sudden_Underflow
                if (k & STRTOG_Inexact) {
                        if (sign)
-                               u->L[_0] = 0x80100000L;
+                               word0(&u[0]) = 0x80100000L;
                        else
-                               u->L[2+_0] = 0x100000L;
+                               word0(&u[1]) = 0x100000L;
                        }
                break;
 #else
@@ -88,80 +84,80 @@ strtodI(CONST char *s, char **sp, double *dd)
 #endif
 
          case STRTOG_Denormal:
-               u->L[_1] = bits[0];
-               u->L[_0] = bits[1];
+               word1(&u[0]) = bits[0];
+               word0(&u[0]) = bits[1];
                goto contain;
 
          case STRTOG_Normal:
-               u->L[_1] = bits[0];
-               u->L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
+               word1(&u[0]) = bits[0];
+               word0(&u[0]) = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
          contain:
                j = k & STRTOG_Inexact;
                if (sign) {
-                       u->L[_0] |= sign;
+                       word0(&u[0]) |= sign;
                        j = STRTOG_Inexact - j;
                        }
                switch(j) {
                  case STRTOG_Inexlo:
 #ifdef Sudden_Underflow
                        if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
-                               u->L[2+_0] = u->L[_0] + 0x3500000;
-                               u->L[2+_1] = u->L[_1];
-                               u->d[1] += ulp(u->d[1]);
-                               u->L[2+_0] -= 0x3500000;
-                               if (!(u->L[2+_0] & 0x7ff00000)) {
-                                       u->L[2+_0] = sign;
-                                       u->L[2+_1] = 0;
+                               word0(&u[1]) = word0(&u[0]) + 0x3500000;
+                               word1(&u[1]) = word1(&u[0]);
+                               dval(&u[1]) += ulp(&u[1]);
+                               word0(&u[1]) -= 0x3500000;
+                               if (!(word0(&u[1]) & 0x7ff00000)) {
+                                       word0(&u[1]) = sign;
+                                       word1(&u[1]) = 0;
                                        }
                                }
                        else
 #endif
-                       u->d[1] = u->d[0] + ulp(u->d[0]);
+                       dval(&u[1]) = dval(&u[0]) + ulp(&u[0]);
                        break;
                  case STRTOG_Inexhi:
-                       u->d[1] = u->d[0];
+                       dval(&u[1]) = dval(&u[0]);
 #ifdef Sudden_Underflow
-                       if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
-                               u->L[_0] += 0x3500000;
-                               u->d[0] -= ulpdown(u->d);
-                               u->L[_0] -= 0x3500000;
-                               if (!(u->L[_0] & 0x7ff00000)) {
-                                       u->L[_0] = sign;
-                                       u->L[_1] = 0;
+                       if ((word0(&u[0]) & 0x7ff00000) < 0x3500000) {
+                               word0(&u[0]) += 0x3500000;
+                               dval(&u[0]) -= ulpdown(u);
+                               word0(&u[0]) -= 0x3500000;
+                               if (!(word0(&u[0]) & 0x7ff00000)) {
+                                       word0(&u[0]) = sign;
+                                       word1(&u[0]) = 0;
                                        }
                                }
                        else
 #endif
-                       u->d[0] -= ulpdown(u->d);
+                       dval(&u[0]) -= ulpdown(u);
                        break;
                  default:
-                       u->d[1] = u->d[0];
+                       dval(&u[1]) = dval(&u[0]);
                  }
                break;
 
          case STRTOG_Infinite:
-               u->L[_0] = u->L[2+_0] = sign | 0x7ff00000;
-               u->L[_1] = u->L[2+_1] = 0;
+               word0(&u[0]) = word0(&u[1]) = sign | 0x7ff00000;
+               word1(&u[0]) = word1(&u[1]) = 0;
                if (k & STRTOG_Inexact) {
                        if (sign) {
-                               u->L[2+_0] = 0xffefffffL;
-                               u->L[2+_1] = 0xffffffffL;
+                               word0(&u[1]) = 0xffefffffL;
+                               word1(&u[1]) = 0xffffffffL;
                                }
                        else {
-                               u->L[_0] = 0x7fefffffL;
-                               u->L[_1] = 0xffffffffL;
+                               word0(&u[0]) = 0x7fefffffL;
+                               word1(&u[0]) = 0xffffffffL;
                                }
                        }
                break;
 
          case STRTOG_NaN:
-               u->L[0] = u->L[2] = d_QNAN0;
-               u->L[1] = u->L[3] = d_QNAN1;
+               u->L[0] = (u+1)->L[0] = d_QNAN0;
+               u->L[1] = (u+1)->L[1] = d_QNAN1;
                break;
 
          case STRTOG_NaNbits:
-               u->L[_0] = u->L[2+_0] = 0x7ff00000 | sign | bits[1];
-               u->L[_1] = u->L[2+_1] = bits[0];
+               word0(&u[0]) = word0(&u[1]) = 0x7ff00000 | sign | bits[1];
+               word1(&u[0]) = word1(&u[1]) = bits[0];
          }
        return k;
        }
index a69ce30..5059869 100644 (file)
@@ -172,9 +172,9 @@ set_ones(Bigint *b, int n)
 rvOK
 #ifdef KR_headers
  (d, fpi, exp, bits, exact, rd, irv)
double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
 #else
- (double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
+ (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
 #endif
 {
        Bigint *b;
@@ -182,7 +182,7 @@ rvOK
        int bdif, e, j, k, k1, nb, rv;
 
        carry = rv = 0;
-       b = d2b(d, &e, &bdif);
+       b = d2b(dval(d), &e, &bdif);
        bdif -= nb = fpi->nbits;
        e += bdif;
        if (bdif <= 0) {
@@ -291,9 +291,9 @@ rvOK
 
  static int
 #ifdef KR_headers
-mantbits(d) double d;
+mantbits(d) U *d;
 #else
-mantbits(double d)
+mantbits(U *d)
 #endif
 {
        ULong L;
@@ -327,8 +327,9 @@ strtodg
        int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
        int sudden_underflow;
        CONST char *s, *s0, *s1;
-       double adj, adj0, rv, tol;
+       double adj0, tol;
        Long L;
+       U adj, rv;
        ULong *b, *be, y, z;
        Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
 #ifdef USE_LOCALE /*{{*/
@@ -341,7 +342,7 @@ strtodg
        static int dplen;
        if (!(s0 = decimalpoint_cache)) {
                s0 = localeconv()->decimal_point;
-               if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
+               if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
                        strcpy(decimalpoint_cache, s0);
                        s0 = decimalpoint_cache;
                        }
@@ -355,7 +356,7 @@ strtodg
 
        irv = STRTOG_Zero;
        denorm = sign = nz0 = nz = 0;
-       dval(rv) = 0.;
+       dval(&rv) = 0.;
        rvb = 0;
        nbits = fpi->nbits;
        for(s = s00;;s++) switch(*s) {
@@ -547,13 +548,13 @@ strtodg
        if (!nd0)
                nd0 = nd;
        k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-       dval(rv) = y;
+       dval(&rv) = y;
        if (k > 9)
-               dval(rv) = tens[k - 9] * dval(rv) + z;
+               dval(&rv) = tens[k - 9] * dval(&rv) + z;
        bd0 = 0;
        if (nbits <= P && nd <= DBL_DIG) {
                if (!e) {
-                       if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
+                       if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
                                goto ret;
                        }
                else if (e > 0) {
@@ -561,9 +562,9 @@ strtodg
 #ifdef VAX
                                goto vax_ovfl_check;
 #else
-                               i = fivesbits[e] + mantbits(dval(rv)) <= P;
-                               /* rv = */ rounded_product(dval(rv), tens[e]);
-                               if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
+                               i = fivesbits[e] + mantbits(&rv) <= P;
+                               /* rv = */ rounded_product(dval(&rv), tens[e]);
+                               if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
                                        goto ret;
                                e1 -= e;
                                goto rv_notOK;
@@ -576,32 +577,32 @@ strtodg
                                 */
                                e2 = e - i;
                                e1 -= i;
-                               dval(rv) *= tens[i];
+                               dval(&rv) *= tens[i];
 #ifdef VAX
                                /* VAX exponent range is so narrow we must
                                 * worry about overflow here...
                                 */
  vax_ovfl_check:
-                               dval(adj) = dval(rv);
-                               word0(adj) -= P*Exp_msk1;
-                               /* adj = */ rounded_product(dval(adj), tens[e2]);
-                               if ((word0(adj) & Exp_mask)
+                               dval(&adj) = dval(&rv);
+                               word0(&adj) -= P*Exp_msk1;
+                               /* adj = */ rounded_product(dval(&adj), tens[e2]);
+                               if ((word0(&adj) & Exp_mask)
                                 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
                                        goto rv_notOK;
-                               word0(adj) += P*Exp_msk1;
-                               dval(rv) = dval(adj);
+                               word0(&adj) += P*Exp_msk1;
+                               dval(&rv) = dval(&adj);
 #else
-                               /* rv = */ rounded_product(dval(rv), tens[e2]);
+                               /* rv = */ rounded_product(dval(&rv), tens[e2]);
 #endif
-                               if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
+                               if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
                                        goto ret;
                                e1 -= e2;
                                }
                        }
 #ifndef Inaccurate_Divide
                else if (e >= -Ten_pmax) {
-                       /* rv = */ rounded_quotient(dval(rv), tens[-e]);
-                       if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
+                       /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
+                       if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
                                goto ret;
                        e1 -= e;
                        }
@@ -615,45 +616,45 @@ strtodg
        e2 = 0;
        if (e1 > 0) {
                if ( (i = e1 & 15) !=0)
-                       dval(rv) *= tens[i];
+                       dval(&rv) *= tens[i];
                if (e1 &= ~15) {
                        e1 >>= 4;
-                       while(e1 >= (1 << n_bigtens-1)) {
-                               e2 += ((word0(rv) & Exp_mask)
+                       while(e1 >= (1 << (n_bigtens-1))) {
+                               e2 += ((word0(&rv) & Exp_mask)
                                        >> Exp_shift1) - Bias;
-                               word0(rv) &= ~Exp_mask;
-                               word0(rv) |= Bias << Exp_shift1;
-                               dval(rv) *= bigtens[n_bigtens-1];
-                               e1 -= 1 << n_bigtens-1;
+                               word0(&rv) &= ~Exp_mask;
+                               word0(&rv) |= Bias << Exp_shift1;
+                               dval(&rv) *= bigtens[n_bigtens-1];
+                               e1 -= 1 << (n_bigtens-1);
                                }
-                       e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
-                       word0(rv) &= ~Exp_mask;
-                       word0(rv) |= Bias << Exp_shift1;
+                       e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
+                       word0(&rv) &= ~Exp_mask;
+                       word0(&rv) |= Bias << Exp_shift1;
                        for(j = 0; e1 > 0; j++, e1 >>= 1)
                                if (e1 & 1)
-                                       dval(rv) *= bigtens[j];
+                                       dval(&rv) *= bigtens[j];
                        }
                }
        else if (e1 < 0) {
                e1 = -e1;
                if ( (i = e1 & 15) !=0)
-                       dval(rv) /= tens[i];
+                       dval(&rv) /= tens[i];
                if (e1 &= ~15) {
                        e1 >>= 4;
-                       while(e1 >= (1 << n_bigtens-1)) {
-                               e2 += ((word0(rv) & Exp_mask)
+                       while(e1 >= (1 << (n_bigtens-1))) {
+                               e2 += ((word0(&rv) & Exp_mask)
                                        >> Exp_shift1) - Bias;
-                               word0(rv) &= ~Exp_mask;
-                               word0(rv) |= Bias << Exp_shift1;
-                               dval(rv) *= tinytens[n_bigtens-1];
-                               e1 -= 1 << n_bigtens-1;
+                               word0(&rv) &= ~Exp_mask;
+                               word0(&rv) |= Bias << Exp_shift1;
+                               dval(&rv) *= tinytens[n_bigtens-1];
+                               e1 -= 1 << (n_bigtens-1);
                                }
-                       e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
-                       word0(rv) &= ~Exp_mask;
-                       word0(rv) |= Bias << Exp_shift1;
+                       e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
+                       word0(&rv) &= ~Exp_mask;
+                       word0(&rv) |= Bias << Exp_shift1;
                        for(j = 0; e1 > 0; j++, e1 >>= 1)
                                if (e1 & 1)
-                                       dval(rv) *= tinytens[j];
+                                       dval(&rv) *= tinytens[j];
                        }
                }
 #ifdef IBM
@@ -664,7 +665,7 @@ strtodg
         */
        e2 <<= 2;
 #endif
-       rvb = d2b(dval(rv), &rve, &rvbits);     /* rv = rvb * 2^rve */
+       rvb = d2b(dval(&rv), &rve, &rvbits);    /* rv = rvb * 2^rve */
        rve += e2;
        if ((j = rvbits - nbits) > 0) {
                rshift(rvb, j);
@@ -842,7 +843,7 @@ strtodg
                                }
                        else
                                irv = STRTOG_Normal | STRTOG_Inexhi;
-                       if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
+                       if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
                                break;
                        if (dsign) {
                                rvb = increment(rvb);
@@ -859,7 +860,7 @@ strtodg
                                }
                        break;
                        }
-               if ((dval(adj) = ratio(delta, bs)) <= 2.) {
+               if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
  adj1:
                        inex = STRTOG_Inexlo;
                        if (dsign) {
@@ -873,15 +874,15 @@ strtodg
                                irv = STRTOG_Underflow | STRTOG_Inexlo;
                                break;
                                }
-                       adj0 = dval(adj) = 1.;
+                       adj0 = dval(&adj) = 1.;
                        }
                else {
-                       adj0 = dval(adj) *= 0.5;
+                       adj0 = dval(&adj) *= 0.5;
                        if (dsign) {
                                asub = 0;
                                inex = STRTOG_Inexlo;
                                }
-                       if (dval(adj) < 2147483647.) {
+                       if (dval(&adj) < 2147483647.) {
                                L = adj0;
                                adj0 -= L;
                                switch(rd) {
@@ -900,12 +901,12 @@ strtodg
                                                inex = STRTOG_Inexact - inex;
                                                }
                                  }
-                               dval(adj) = L;
+                               dval(&adj) = L;
                                }
                        }
                y = rve + rvbits;
 
-               /* adj *= ulp(dval(rv)); */
+               /* adj *= ulp(dval(&rv)); */
                /* if (asub) rv -= adj; else rv += adj; */
 
                if (!denorm && rvbits < nbits) {
@@ -913,7 +914,7 @@ strtodg
                        rve -= j;
                        rvbits = nbits;
                        }
-               ab = d2b(dval(adj), &abe, &abits);
+               ab = d2b(dval(&adj), &abe, &abits);
                if (abe < 0)
                        rshift(ab, -abe);
                else if (abe > 0)
@@ -967,15 +968,15 @@ strtodg
                z = rve + rvbits;
                if (y == z && L) {
                        /* Can we stop now? */
-                       tol = dval(adj) * 5e-16; /* > max rel error */
-                       dval(adj) = adj0 - .5;
-                       if (dval(adj) < -tol) {
+                       tol = dval(&adj) * 5e-16; /* > max rel error */
+                       dval(&adj) = adj0 - .5;
+                       if (dval(&adj) < -tol) {
                                if (adj0 > tol) {
                                        irv |= inex;
                                        break;
                                        }
                                }
-                       else if (dval(adj) > tol && adj0 < 1. - tol) {
+                       else if (dval(&adj) > tol && adj0 < 1. - tol) {
                                irv |= inex;
                                break;
                                }
index f8111b7..a8beb35 100644 (file)
@@ -58,7 +58,7 @@ strtof(CONST char *s, char **sp)
 
          case STRTOG_Normal:
          case STRTOG_NaNbits:
-               u.L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+               u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
                break;
 
          case STRTOG_Denormal:
index c665976..738372d 100644 (file)
@@ -67,8 +67,8 @@ strtopdd(CONST char *s, char **sp, double *dd)
 
          case STRTOG_Normal:
                u->L[_1] = (bits[1] >> 21 | bits[2] << 11) & 0xffffffffL;
-               u->L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
-                         | exp + 0x3ff + 105 << 20;
+               u->L[_0] = (bits[2] >> 21) | ((bits[3] << 11) & 0xfffff)
+                         | ((exp + 0x3ff + 105) << 20);
                exp += 0x3ff + 52;
                if (bits[1] &= 0x1fffff) {
                        i = hi0bits(bits[1]) - 11;
@@ -79,7 +79,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
                        else
                                exp -= i;
                        if (i > 0) {
-                               bits[1] = bits[1] << i | bits[0] >> 32-i;
+                               bits[1] = bits[1] << i | bits[0] >> (32-i);
                                bits[0] = bits[0] << i & 0xffffffffL;
                                }
                        }
@@ -92,11 +92,11 @@ strtopdd(CONST char *s, char **sp, double *dd)
                        else
                                exp -= i;
                        if (i < 32) {
-                               bits[1] = bits[0] >> 32 - i;
+                               bits[1] = bits[0] >> (32 - i);
                                bits[0] = bits[0] << i & 0xffffffffL;
                                }
                        else {
-                               bits[1] = bits[0] << i - 32;
+                               bits[1] = bits[0] << (i - 32);
                                bits[0] = 0;
                                }
                        }
@@ -105,7 +105,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
                        break;
                        }
                u->L[2+_1] = bits[0];
-               u->L[2+_0] = bits[1] & 0xfffff | exp << 20;
+               u->L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
                break;
 
          case STRTOG_Denormal:
@@ -124,10 +124,10 @@ strtopdd(CONST char *s, char **sp, double *dd)
          nearly_normal:
                i = hi0bits(bits[3]) - 11;      /* i >= 12 */
                j = 32 - i;
-               u->L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
-                       | 65 - i << 20;
+               u->L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
+                       | ((65 - i) << 20);
                u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-               u->L[2+_0] = bits[1] & (1L << j) - 1;
+               u->L[2+_0] = bits[1] & ((1L << j) - 1);
                u->L[2+_1] = bits[0];
                break;
 
@@ -136,34 +136,34 @@ strtopdd(CONST char *s, char **sp, double *dd)
                if (i < 0) {
                        j = -i;
                        i += 32;
-                       u->L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
-                       u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-                       u->L[2+_0] = bits[1] & (1L << j) - 1;
+                       u->L[_0] = (bits[2] >> j & 0xfffff) | (33 + j) << 20;
+                       u->L[_1] = ((bits[2] << i) | (bits[1] >> j)) & 0xffffffffL;
+                       u->L[2+_0] = bits[1] & ((1L << j) - 1);
                        u->L[2+_1] = bits[0];
                        break;
                        }
                if (i == 0) {
-                       u->L[_0] = bits[2] & 0xfffff | 33 << 20;
+                       u->L[_0] = (bits[2] & 0xfffff) | (33 << 20);
                        u->L[_1] = bits[1];
                        u->L[2+_0] = 0;
                        u->L[2+_1] = bits[0];
                        break;
                        }
                j = 32 - i;
-               u->L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
-                               | j + 1 << 20;
+               u->L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
+                               | ((j + 1) << 20);
                u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
                u->L[2+_0] = 0;
-               u->L[2+_1] = bits[0] & (1L << j) - 1;
+               u->L[2+_1] = bits[0] & ((1L << j) - 1);
                break;
 
          hardly_normal:
                j = 11 - hi0bits(bits[1]);
                i = 32 - j;
-               u->L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
+               u->L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
                u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
                u->L[2+_0] = 0;
-               u->L[2+_1] = bits[0] & (1L << j) - 1;
+               u->L[2+_1] = bits[0] & ((1L << j) - 1);
                break;
 
          case STRTOG_Infinite:
index 74694a0..23ca5cb 100644 (file)
@@ -58,7 +58,7 @@ strtopf(CONST char *s, char **sp, float *f)
 
          case STRTOG_Normal:
          case STRTOG_NaNbits:
-               L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+               L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
                break;
 
          case STRTOG_Denormal:
index 07d0458..f7a25ff 100644 (file)
@@ -92,7 +92,8 @@ strtopx(CONST char *s, char **sp, void *V)
 
          case STRTOG_Infinite:
                L[_0] = 0x7fff;
-               L[_1] = L[_2] = L[_3] = L[_4] = 0;
+               L[_1] = 0x8000;
+               L[_2] = L[_3] = L[_4] = 0;
                break;
 
          case STRTOG_NaN:
index b0f5cdd..6519a41 100644 (file)
@@ -82,7 +82,8 @@ strtopxL(CONST char *s, char **sp, void *V)
 
          case STRTOG_Infinite:
                L[_0] = 0x7fff << 16;
-               L[_1] = L[_2] = 0;
+               L[_1] = 0x80000000;
+               L[_2] = 0;
                break;
 
          case STRTOG_NaN:
index 9c2b46e..e0b8e6a 100644 (file)
@@ -48,8 +48,8 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 
          case STRTOG_Normal:
                L[_1] = (bits[1] >> 21 | bits[2] << 11) & (ULong)0xffffffffL;
-               L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
-                         | exp + 0x3ff + 105 << 20;
+               L[_0] = (bits[2] >> 21) | (bits[3] << 11 & 0xfffff)
+                         | ((exp + 0x3ff + 105) << 20);
                exp += 0x3ff + 52;
                if (bits[1] &= 0x1fffff) {
                        i = hi0bits(bits[1]) - 11;
@@ -60,7 +60,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
                        else
                                exp -= i;
                        if (i > 0) {
-                               bits[1] = bits[1] << i | bits[0] >> 32-i;
+                               bits[1] = bits[1] << i | bits[0] >> (32-i);
                                bits[0] = bits[0] << i & (ULong)0xffffffffL;
                                }
                        }
@@ -73,11 +73,11 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
                        else
                                exp -= i;
                        if (i < 32) {
-                               bits[1] = bits[0] >> 32 - i;
+                               bits[1] = bits[0] >> (32 - i);
                                bits[0] = bits[0] << i & (ULong)0xffffffffL;
                                }
                        else {
-                               bits[1] = bits[0] << i - 32;
+                               bits[1] = bits[0] << (i - 32);
                                bits[0] = 0;
                                }
                        }
@@ -86,7 +86,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
                        break;
                        }
                L[2+_1] = bits[0];
-               L[2+_0] = bits[1] & 0xfffff | exp << 20;
+               L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
                break;
 
          case STRTOG_Denormal:
@@ -105,10 +105,10 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
          nearly_normal:
                i = hi0bits(bits[3]) - 11;      /* i >= 12 */
                j = 32 - i;
-               L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
-                       | 65 - i << 20;
+               L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
+                       | ((65 - i) << 20);
                L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-               L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
+               L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
                L[2+_1] = bits[0];
                break;
 
@@ -117,34 +117,34 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
                if (i < 0) {
                        j = -i;
                        i += 32;
-                       L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
+                       L[_0] = (bits[2] >> j & 0xfffff) | ((33 + j) << 20);
                        L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-                       L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
+                       L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
                        L[2+_1] = bits[0];
                        break;
                        }
                if (i == 0) {
-                       L[_0] = bits[2] & 0xfffff | 33 << 20;
+                       L[_0] = (bits[2] & 0xfffff) | (33 << 20);
                        L[_1] = bits[1];
                        L[2+_0] = 0;
                        L[2+_1] = bits[0];
                        break;
                        }
                j = 32 - i;
-               L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
-                               | j + 1 << 20;
+               L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
+                               | ((j + 1) << 20);
                L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
                L[2+_0] = 0;
-               L[2+_1] = bits[0] & (1L << j) - 1;
+               L[2+_1] = bits[0] & ((1L << j) - 1);
                break;
 
          hardly_normal:
                j = 11 - hi0bits(bits[1]);
                i = 32 - j;
-               L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
+               L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
                L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
                L[2+_0] = 0;
-               L[2+_1] = bits[0] & ((ULong)1L << j) - 1;
+               L[2+_1] = bits[0] & (((ULong)1L << j) - 1);
                break;
 
          case STRTOG_Infinite:
index 46b0ba2..21c6d28 100644 (file)
@@ -46,7 +46,7 @@ ULtof(ULong *L, ULong *bits, Long exp, int k)
 
          case STRTOG_Normal:
          case STRTOG_NaNbits:
-               L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+               L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
                break;
 
          case STRTOG_Denormal:
index 23f721a..d2f9472 100644 (file)
@@ -80,7 +80,8 @@ ULtox(UShort *L, ULong *bits, Long exp, int k)
 
          case STRTOG_Infinite:
                L[_0] = 0x7fff;
-               L[_1] = L[_2] = L[_3] = L[_4] = 0;
+               L[_1] = 0x8000;
+               L[_2] = L[_3] = L[_4] = 0;
                break;
 
          case STRTOG_NaN:
index ff62a61..58acc80 100644 (file)
@@ -70,7 +70,8 @@ ULtoxL(ULong *L, ULong *bits, Long exp, int k)
 
          case STRTOG_Infinite:
                L[_0] = 0x7fff << 16;
-               L[_1] = L[_2] = 0;
+               L[_1] = 0x80000000;
+               L[_2] = 0;
                break;
 
          case STRTOG_NaN:
index 7810a5c..17e9f86 100644 (file)
@@ -34,13 +34,13 @@ THIS SOFTWARE.
  double
 ulp
 #ifdef KR_headers
-       (x) double x;
+       (x) U *x;
 #else
-       (double x)
+       (U *x)
 #endif
 {
        Long L;
-       double a;
+       U a;
 
        L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
 #ifndef Sudden_Underflow
@@ -49,22 +49,22 @@ ulp
 #ifdef IBM
                L |= Exp_msk1 >> 4;
 #endif
-               word0(a) = L;
-               word1(a) = 0;
+               word0(&a) = L;
+               word1(&a) = 0;
 #ifndef Sudden_Underflow
                }
        else {
                L = -L >> Exp_shift;
                if (L < Exp_shift) {
-                       word0(a) = 0x80000 >> L;
-                       word1(a) = 0;
+                       word0(&a) = 0x80000 >> L;
+                       word1(&a) = 0;
                        }
                else {
-                       word0(a) = 0;
+                       word0(&a) = 0;
                        L -= Exp_shift;
-                       word1(a) = L >= 31 ? 1 : 1 << 31 - L;
+                       word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
                        }
                }
 #endif
-       return a;
+       return dval(&a);
        }