From cd3f5db16a8866d7ed2e97a05ab74e1f631feaf2 Mon Sep 17 00:00:00 2001 From: Braden Obrzut Date: Sat, 4 Jan 2014 16:58:21 -0500 Subject: [PATCH] - Updated gdtoa in order to fix clang warnings. --- gdtoa/README | 84 ++++++- gdtoa/arithchk.c | 2 +- gdtoa/dmisc.c | 2 +- gdtoa/dtoa.c | 210 +++++++++-------- gdtoa/g_Qfmt.c | 17 +- gdtoa/g__fmt.c | 126 +++++++++- gdtoa/g_ddfmt.c | 47 ++-- gdtoa/g_dfmt.c | 26 +- gdtoa/g_ffmt.c | 17 +- gdtoa/g_xLfmt.c | 17 +- gdtoa/g_xfmt.c | 23 +- gdtoa/gdtoa.c | 196 +++++++-------- gdtoa/gdtoa.h | 42 ++-- gdtoa/gdtoa_fltrnds.h | 18 ++ gdtoa/gdtoaimp.h | 81 +++++-- gdtoa/gethex.c | 155 +++++++++--- gdtoa/gmisc.c | 2 +- gdtoa/hd_init.c | 24 +- gdtoa/hexnan.c | 27 ++- gdtoa/misc.c | 136 ++++------- gdtoa/qnan.c | 45 ++-- gdtoa/smisc.c | 26 +- gdtoa/strtoIQ.c | 2 +- gdtoa/strtoId.c | 2 +- gdtoa/strtoIdd.c | 4 +- gdtoa/strtoIf.c | 2 +- gdtoa/strtoIg.c | 20 +- gdtoa/strtoIx.c | 2 +- gdtoa/strtoIxL.c | 2 +- gdtoa/strtod.c | 536 ++++++++++++++++++++++++------------------ gdtoa/strtodI.c | 90 ++++--- gdtoa/strtodg.c | 203 ++++++++++------ gdtoa/strtodnrp.c | 2 +- gdtoa/strtof.c | 11 +- gdtoa/strtopQ.c | 20 +- gdtoa/strtopd.c | 9 +- gdtoa/strtopdd.c | 47 ++-- gdtoa/strtopf.c | 11 +- gdtoa/strtopx.c | 24 +- gdtoa/strtopxL.c | 20 +- gdtoa/strtorQ.c | 15 +- gdtoa/strtord.c | 11 +- gdtoa/strtordd.c | 58 ++--- gdtoa/strtorf.c | 11 +- gdtoa/strtorx.c | 20 +- gdtoa/strtorxL.c | 18 +- gdtoa/ulp.c | 21 +- 47 files changed, 1528 insertions(+), 956 deletions(-) create mode 100644 gdtoa/gdtoa_fltrnds.h diff --git a/gdtoa/README b/gdtoa/README index cf1947aa2..1bf7d91e4 100644 --- a/gdtoa/README +++ b/gdtoa/README @@ -56,7 +56,9 @@ two letters: whose sum is the desired value For decimal -> binary conversions, there are three families of -helper routines: one for round-nearest: +helper routines: one for round-nearest (or the current rounding +mode on IEEE-arithmetic systems that provide the C99 fegetround() +function, if compiled with -DHonor_FLT_ROUNDS): strtof strtod @@ -150,12 +152,14 @@ suffer double rounding due to use of extended-precision registers. For some conversions this variant of strtod is less efficient than the one in strtod.c when the latter is run with 53-bit rounding precision. -The values that the strto* routines return for NaNs are determined by -gd_qnan.h, which the makefile generates by running the program whose -source is qnan.c. Note that the rules for distinguishing signaling -from quiet NaNs are system-dependent. For cross-compilation, you need -to determine arith.h and gd_qnan.h suitably, e.g., using the -arithmetic of the target machine. +When float or double are involved, the values that the strto* routines +return for NaNs are determined by gd_qnan.h, which the makefile +generates by running the program whose source is qnan.c. For other +types, default NaN values are specified in g__fmt.c and may need +adjusting. Note that the rules for distinguishing signaling from +quiet NaNs are system-dependent. For cross-compilation, you need to +determine arith.h and gd_qnan.h suitably, e.g., using the arithmetic +of the target machine. C99's hexadecimal floating-point constants are recognized by the strto* routines (but this feature has not yet been heavily tested). @@ -170,10 +174,11 @@ hexadecimal digits, it is taken for the fraction bits of the resulting NaN; if there are two or more strings of hexadecimal digits, each string is assigned to the next available sequence of 32-bit words of fractions bits (starting with the most significant), right-aligned in -each sequence. +each sequence. Strings of hexadecimal digits may be preceded by "0x" +or "0X". -For binary -> decimal conversions, I've provided just one family -of helper routines: +For binary -> decimal conversions, I've provided a family of helper +routines: g_ffmt g_dfmt @@ -181,6 +186,12 @@ of helper routines: g_xfmt g_xLfmt g_Qfmt + g_ffmt_p + g_dfmt_p + g_ddfmt_p + g_xfmt_p + g_xLfmt_p + g_Qfmt_p which do a "%g" style conversion either to a specified number of decimal places (if their ndig argument is positive), or to the shortest @@ -191,6 +202,36 @@ in the buffer, if the buffer was long enough, or 0. Other forms of conversion are easily done with the help of gdtoa(), such as %e or %f style and conversions with direction of rounding specified (so that, if desired, the decimal value is either >= or <= the binary value). +On IEEE-arithmetic systems that provide the C99 fegetround() function, +if compiled with -DHonor_FLT_ROUNDS, these routines honor the current +rounding mode. For pedants, the ...fmt_p() routines are similar to the +...fmt() routines, but have an additional final int argument, nik, +that for conversions of Infinity or NaN, determines whether upper, +lower, or mixed case is used, whether (...) is added to NaN values, +and whether the sign of a NaN is reported or suppressed: + + nik = ic + 6*(nb + 3*ns), + +where ic with 0 <= ic < 6 controls the rendering of Infinity and NaN: + + 0 ==> Infinity or NaN + 1 ==> infinity or nan + 2 ==> INFINITY or NAN + 3 ==> Inf or NaN + 4 ==> inf or nan + 5 ==> INF or NAN + +nb with 0 <= nb < 3 determines whether NaN values are rendered +as NaN(...): + + 0 ==> no + 1 ==> yes + 2 ==> no for default NaN values; yes otherwise + +ns = 0 or 1 determines whether the sign of NaN values reported: + + 0 ==> distinguish NaN and -NaN + 1 ==> report both as NaN For an example of more general conversions based on dtoa(), see netlib's "printf.c from ampl/solvers". @@ -332,5 +373,28 @@ Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes the decimal-point character to be taken from the current locale; otherwise it is '.'. +Source files dtoa.c and strtod.c in this directory are derived from +netlib's "dtoa.c from fp" and are meant to function equivalently. +When compiled with Honor_FLT_ROUNDS #defined (on systems that provide +FLT_ROUNDS and fegetround() as specified in the C99 standard), they +honor the current rounding mode. Because FLT_ROUNDS is buggy on some +(Linux) systems -- not reflecting calls on fesetround(), as the C99 +standard says it should -- when Honor_FLT_ROUNDS is #defined, the +current rounding mode is obtained from fegetround() rather than from +FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined. + +Compile with -DUSE_LOCALE to use the current locale; otherwise +decimal points are assumed to be '.'. With -DUSE_LOCALE, unless +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 "."). diff --git a/gdtoa/arithchk.c b/gdtoa/arithchk.c index 3211aeda4..ef6cda3db 100644 --- a/gdtoa/arithchk.c +++ b/gdtoa/arithchk.c @@ -106,7 +106,7 @@ ccheck() long Cray1; /* Cray1 = 4617762693716115456 -- without overflow on non-Crays */ - Cray1 = printf(emptyfmt) < 0 ? 0 : 4617762; + Cray1 = printf("%s", emptyfmt) < 0 ? 0 : 4617762; if (printf(emptyfmt, Cray1) >= 0) Cray1 = 1000000*Cray1 + 693716; if (printf(emptyfmt, Cray1) >= 0) diff --git a/gdtoa/dmisc.c b/gdtoa/dmisc.c index 72d271ca9..3e712511b 100644 --- a/gdtoa/dmisc.c +++ b/gdtoa/dmisc.c @@ -46,7 +46,7 @@ rv_alloc(int i) j = sizeof(ULong); for(k = 0; - sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i; + sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)(i); j <<= 1) k++; r = (int*)Balloc(k); diff --git a/gdtoa/dtoa.c b/gdtoa/dtoa.c index 4fc6633b7..c96e6a545 100644 --- a/gdtoa/dtoa.c +++ b/gdtoa/dtoa.c @@ -30,7 +30,6 @@ THIS SOFTWARE. * with " at " changed at "@" and " dot " changed to "."). */ #include "gdtoaimp.h" -#include /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. * @@ -67,7 +66,6 @@ THIS SOFTWARE. */ #ifdef Honor_FLT_ROUNDS -#define Rounding rounding #undef Check_FLT_ROUNDS #define Check_FLT_ROUNDS #else @@ -77,17 +75,17 @@ 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 of ecvt and fcvt; trailing zeros are suppressed from the returned string. If not null, *rve is set to point to the end of the return value. If d is +-Infinity or NaN, - then *decpt is set to INT_MAX. + then *decpt is set to 9999. mode: 0 ==> shortest string that yields d when read in @@ -129,12 +127,22 @@ dtoa U d, d2, eps; double ds; char *s, *s0; -#ifdef Honor_FLT_ROUNDS - int rounding; -#endif #ifdef SET_INEXACT int inexact, oldinexact; #endif +#ifdef Honor_FLT_ROUNDS /*{*/ + int Rounding; +#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ + Rounding = Flt_Rounds; +#else /*}{*/ + Rounding = 1; + switch(fegetround()) { + case FE_TOWARDZERO: Rounding = 0; break; + case FE_UPWARD: Rounding = 2; break; + case FE_DOWNWARD: Rounding = 3; + } +#endif /*}}*/ +#endif /*}*/ #ifndef MULTIPLE_THREADS if (dtoa_result) { @@ -142,36 +150,35 @@ dtoa dtoa_result = 0; } #endif - - dval(d) = _d; - 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 = INT_MAX; + *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); } @@ -181,35 +188,35 @@ dtoa inexact = 1; #endif #ifdef Honor_FLT_ROUNDS - if ((rounding = Flt_Rounds) >= 2) { + if (Rounding >= 2) { if (*sign) - rounding = rounding == 2 ? 0 : 2; + Rounding = Rounding == 2 ? 0 : 2; else - if (rounding != 2) - rounding = 0; + if (Rounding != 2) + Rounding = 0; } #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 ); @@ -238,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; } @@ -291,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; @@ -319,7 +327,7 @@ dtoa s = s0 = rv_alloc(i); #ifdef Honor_FLT_ROUNDS - if (mode > 1 && rounding != 1) + if (mode > 1 && Rounding != 1) leftright = 0; #endif @@ -328,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 */ @@ -338,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++) @@ -346,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; } @@ -380,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 = (Long)dval(d); - dval(d) -= L; + L = (Long)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; @@ -420,7 +428,7 @@ dtoa #endif fast_failed: s = s0; - dval(d) = dval(d2); + dval(&d) = dval(&d2); k = k0; ilim = ilim0; } @@ -432,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 @@ -456,13 +464,18 @@ dtoa if (i == ilim) { #ifdef Honor_FLT_ROUNDS if (mode > 1) - switch(rounding) { + switch(Rounding) { case 0: goto ret1; 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) { @@ -524,12 +537,12 @@ dtoa spec_case = 0; if ((mode < 2 || leftright) #ifdef Honor_FLT_ROUNDS - && rounding == 1 + && 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 */ @@ -615,9 +628,9 @@ 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 + && Rounding >= 1 #endif ) { if (dig == '9') @@ -632,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; @@ -645,7 +658,7 @@ dtoa } #ifdef Honor_FLT_ROUNDS if (mode > 1) - switch(rounding) { + switch(Rounding) { case 0: goto accept_dig; case 2: goto keep_dig; } @@ -653,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; } @@ -663,7 +680,7 @@ dtoa } if (j1 > 0) { #ifdef Honor_FLT_ROUNDS - if (!rounding) + if (!Rounding) goto accept_dig; #endif if (dig == '9') { /* possible if i == 1 */ @@ -706,14 +723,19 @@ dtoa /* Round off last digit */ #ifdef Honor_FLT_ROUNDS - switch(rounding) { + switch(Rounding) { case 0: goto trimzeros; case 2: goto roundoff; } #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) { @@ -724,7 +746,9 @@ dtoa ++*s++; } else { +#ifdef Honor_FLT_ROUNDS trimzeros: +#endif while(*--s == '0'); s++; } @@ -739,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) diff --git a/gdtoa/g_Qfmt.c b/gdtoa/g_Qfmt.c index 9d227d170..0f0697005 100644 --- a/gdtoa/g_Qfmt.c +++ b/gdtoa/g_Qfmt.c @@ -51,19 +51,24 @@ THIS SOFTWARE. char* #ifdef KR_headers -g_Qfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize; +g_Qfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize; #else -g_Qfmt(char *buf, void *V, int ndig, unsigned bufsize) +g_Qfmt(char *buf, void *V, int ndig, size_t bufsize) #endif { - static CONST FPI fpi = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0 }; + static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0, Int_max }; char *b, *s, *se; ULong bits[4], *L, sign; int decpt, ex, i, mode; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif if (ndig < 0) ndig = 0; - if (bufsize < (unsigned)(ndig + 10)) + if (bufsize < (size_t)(ndig + 10)) return 0; L = (ULong*)V; @@ -109,6 +114,6 @@ g_Qfmt(char *buf, void *V, int ndig, unsigned bufsize) return 0; mode = 0; } - s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - return g__fmt(buf, s, se, decpt, sign); + s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se); + return g__fmt(buf, s, se, decpt, sign, bufsize); } diff --git a/gdtoa/g__fmt.c b/gdtoa/g__fmt.c index 021ecfb57..652c82b68 100644 --- a/gdtoa/g__fmt.c +++ b/gdtoa/g__fmt.c @@ -35,26 +35,77 @@ THIS SOFTWARE. #include "locale.h" #endif +#ifndef ldus_QNAN0 +#define ldus_QNAN0 0x7fff +#endif +#ifndef ldus_QNAN1 +#define ldus_QNAN1 0xc000 +#endif +#ifndef ldus_QNAN2 +#define ldus_QNAN2 0 +#endif +#ifndef ldus_QNAN3 +#define ldus_QNAN3 0 +#endif +#ifndef ldus_QNAN4 +#define ldus_QNAN4 0 +#endif + + const char *InfName[6] = { "Infinity", "infinity", "INFINITY", "Inf", "inf", "INF" }; + const char *NanName[3] = { "NaN", "nan", "NAN" }; + ULong NanDflt_Q_D2A[4] = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff }; + ULong NanDflt_d_D2A[2] = { d_QNAN1, d_QNAN0 }; + ULong NanDflt_f_D2A[1] = { f_QNAN }; + ULong NanDflt_xL_D2A[3] = { 1, 0x80000000, 0x7fff0000 }; + UShort NanDflt_ldus_D2A[5] = { ldus_QNAN4, ldus_QNAN3, ldus_QNAN2, ldus_QNAN1, ldus_QNAN0 }; + char * #ifdef KR_headers -g__fmt(b, s, se, decpt, sign) char *b; char *s; char *se; int decpt; ULong sign; +g__fmt(b, s, se, decpt, sign, blen) char *b; char *s; char *se; int decpt; ULong sign; size_t blen; #else -g__fmt(char *b, char *s, char *se, int decpt, ULong sign) +g__fmt(char *b, char *s, char *se, int decpt, ULong sign, size_t blen) #endif { int i, j, k; - char *s0 = s; + char *be, *s0; + size_t len; #ifdef USE_LOCALE - char decimalpoint = *localeconv()->decimal_point; +#ifdef NO_LOCALE_CACHE + char *decimalpoint = localeconv()->decimal_point; + size_t dlen = strlen(decimalpoint); #else -#define decimalpoint '.' + char *decimalpoint; + static char *decimalpoint_cache; + static size_t dlen; + if (!(s0 = decimalpoint_cache)) { + s0 = localeconv()->decimal_point; + dlen = strlen(s0); + if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { + strcpy(decimalpoint_cache, s0); + s0 = decimalpoint_cache; + } + } + decimalpoint = s0; #endif +#else +#define dlen 0 +#endif + s0 = s; + len = (se-s) + dlen + 6; /* 6 = sign + e+dd + trailing null */ + if (blen < len) + goto ret0; + be = b + blen - 1; if (sign) *b++ = '-'; if (decpt <= -4 || decpt > se - s + 5) { *b++ = *s++; if (*s) { - *b++ = decimalpoint; +#ifdef USE_LOCALE + while((*b = *decimalpoint++)) + ++b; +#else + *b++ = '.'; +#endif while((*b = *s++) !=0) b++; } @@ -69,6 +120,8 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign) for(j = 2, k = 10; 10*k <= decpt; j++, k *= 10){} for(;;) { i = decpt / k; + if (b >= be) + goto ret0; *b++ = i + '0'; if (--j <= 0) break; @@ -78,22 +131,73 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign) *b = 0; } else if (decpt <= 0) { - *b++ = decimalpoint; +#ifdef USE_LOCALE + while((*b = *decimalpoint++)) + ++b; +#else + *b++ = '.'; +#endif + if (be < b - decpt + (se - s)) + goto ret0; for(; decpt < 0; decpt++) *b++ = '0'; - while((*b = *s++) !=0) + while((*b = *s++) != 0) b++; } else { - while((*b = *s++) !=0) { + while((*b = *s++) != 0) { b++; - if (--decpt == 0 && *s) - *b++ = decimalpoint; + if (--decpt == 0 && *s) { +#ifdef USE_LOCALE + while(*b = *decimalpoint++) + ++b; +#else + *b++ = '.'; +#endif + } + } + if (b + decpt > be) { + ret0: + b = 0; + goto ret; } for(; decpt > 0; decpt--) *b++ = '0'; *b = 0; } + ret: freedtoa(s0); return b; } + + char * +add_nanbits_D2A(char *b, size_t blen, ULong *bits, int nb) +{ + ULong t; + char *rv; + int i, j; + size_t L; + static char Hexdig[16] = "0123456789abcdef"; + + while(!bits[--nb]) + if (!nb) + return b; + L = 8*nb + 3; + t = bits[nb]; + do ++L; while((t >>= 4)); + if (L > blen) + return b; + b += L; + *--b = 0; + rv = b; + *--b = /*(*/ ')'; + for(i = 0; i < nb; ++i) { + t = bits[i]; + for(j = 0; j < 8; ++j, t >>= 4) + *--b = Hexdig[t & 0xf]; + } + t = bits[nb]; + do *--b = Hexdig[t & 0xf]; while(t >>= 4); + *--b = '('; /*)*/ + return rv; + } diff --git a/gdtoa/g_ddfmt.c b/gdtoa/g_ddfmt.c index d96e78dbe..5ce4a076b 100644 --- a/gdtoa/g_ddfmt.c +++ b/gdtoa/g_ddfmt.c @@ -33,9 +33,9 @@ THIS SOFTWARE. char * #ifdef KR_headers -g_ddfmt(buf, dd, ndig, bufsize) char *buf; double *dd; int ndig; unsigned 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, unsigned bufsize) +g_ddfmt(char *buf, double *dd0, int ndig, size_t bufsize) #endif { FPI fpi; @@ -43,12 +43,28 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned 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! */ + Rounding = Flt_Rounds; +#else /*}{*/ + Rounding = 1; + switch(fegetround()) { + case FE_TOWARDZERO: Rounding = 0; break; + case FE_UPWARD: Rounding = 2; break; + case FE_DOWNWARD: Rounding = 3; + } +#endif /*}}*/ +#else /*}{*/ +#define Rounding FPI_Round_near +#endif /*}}*/ - if (bufsize < 10 || bufsize < (unsigned)(ndig + 8)) + if (bufsize < 10 || bufsize < (size_t)(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]) { @@ -73,7 +89,7 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned 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) @@ -84,16 +100,16 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned 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); @@ -136,7 +152,7 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize) } mode = 2; if (ndig <= 0) { - if (bufsize < (unsigned)((int)(fpi.nbits * .301029995664) + 10)) { + if (bufsize < (size_t)((int)(fpi.nbits * .301029995664) + 10)) { Bfree(z); return 0; } @@ -144,11 +160,12 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize) } fpi.emin = 1-1023-53+1; fpi.emax = 2046-1023-106+1; - fpi.rounding = FPI_Round_near; + fpi.rounding = Rounding; fpi.sudden_underflow = 0; + fpi.int_max = Int_max; i = STRTOG_Normal; s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - b = g__fmt(buf, s, se, decpt, z->sign); + b = g__fmt(buf, s, se, decpt, z->sign, bufsize); Bfree(z); return b; } diff --git a/gdtoa/g_dfmt.c b/gdtoa/g_dfmt.c index 74b4ccc28..d8e1438c4 100644 --- a/gdtoa/g_dfmt.c +++ b/gdtoa/g_dfmt.c @@ -33,25 +33,32 @@ THIS SOFTWARE. char* #ifdef KR_headers -g_dfmt(buf, d, ndig, bufsize) char *buf; double *d; int ndig; unsigned bufsize; +g_dfmt(buf, d, ndig, bufsize) char *buf; double *d; int ndig; size_t bufsize; #else -g_dfmt(char *buf, double *d, int ndig, unsigned bufsize) +g_dfmt(char *buf, double *d, int ndig, size_t bufsize) #endif { - static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0 }; + static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0, Int_max }; char *b, *s, *se; ULong bits[2], *L, sign; int decpt, ex, i, mode; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif if (ndig < 0) ndig = 0; - if (bufsize < (unsigned)(ndig + 10)) + if (bufsize < (size_t)(ndig + 10)) return 0; L = (ULong*)d; sign = L[_0] & 0x80000000L; if ((L[_0] & 0x7ff00000) == 0x7ff00000) { /* Infinity or NaN */ + if (bufsize < 10) + return 0; if (L[_0] & 0xfffff || L[_1]) { return strcp(buf, "NaN"); } @@ -78,12 +85,11 @@ g_dfmt(char *buf, double *d, int ndig, unsigned bufsize) ex = 1; ex -= 0x3ff + 52; mode = 2; - if (ndig <= 0) { - if (bufsize < 25) - return 0; + if (ndig <= 0) mode = 0; - } i = STRTOG_Normal; - s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - return g__fmt(buf, s, se, decpt, sign); + 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); } diff --git a/gdtoa/g_ffmt.c b/gdtoa/g_ffmt.c index 9b2c5fea6..30b53ae7e 100644 --- a/gdtoa/g_ffmt.c +++ b/gdtoa/g_ffmt.c @@ -33,19 +33,24 @@ THIS SOFTWARE. char* #ifdef KR_headers -g_ffmt(buf, f, ndig, bufsize) char *buf; float *f; int ndig; unsigned bufsize; +g_ffmt(buf, f, ndig, bufsize) char *buf; float *f; int ndig; size_t bufsize; #else -g_ffmt(char *buf, float *f, int ndig, unsigned bufsize) +g_ffmt(char *buf, float *f, int ndig, size_t bufsize) #endif { - static CONST FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, 0 }; + static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, 0, 6 }; char *b, *s, *se; ULong bits[1], *L, sign; int decpt, ex, i, mode; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif if (ndig < 0) ndig = 0; - if (bufsize < (unsigned)(ndig + 10)) + if (bufsize < (size_t)(ndig + 10)) return 0; L = (ULong*)f; @@ -83,6 +88,6 @@ g_ffmt(char *buf, float *f, int ndig, unsigned bufsize) mode = 0; } i = STRTOG_Normal; - s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - return g__fmt(buf, s, se, decpt, sign); + s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se); + return g__fmt(buf, s, se, decpt, sign, bufsize); } diff --git a/gdtoa/g_xLfmt.c b/gdtoa/g_xLfmt.c index a2a993b59..5cda8d59e 100644 --- a/gdtoa/g_xLfmt.c +++ b/gdtoa/g_xLfmt.c @@ -49,19 +49,24 @@ THIS SOFTWARE. char* #ifdef KR_headers -g_xLfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize; +g_xLfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize; #else -g_xLfmt(char *buf, void *V, int ndig, unsigned bufsize) +g_xLfmt(char *buf, void *V, int ndig, size_t bufsize) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 }; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0, Int_max }; char *b, *s, *se; ULong bits[2], *L, sign; int decpt, ex, i, mode; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif if (ndig < 0) ndig = 0; - if (bufsize < (unsigned)(ndig + 10)) + if (bufsize < (size_t)(ndig + 10)) return 0; L = (ULong*)V; @@ -103,6 +108,6 @@ g_xLfmt(char *buf, void *V, int ndig, unsigned bufsize) return 0; mode = 0; } - s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - return g__fmt(buf, s, se, decpt, sign); + s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se); + return g__fmt(buf, s, se, decpt, sign, bufsize); } diff --git a/gdtoa/g_xfmt.c b/gdtoa/g_xfmt.c index c51a16dce..a0baa518c 100644 --- a/gdtoa/g_xfmt.c +++ b/gdtoa/g_xfmt.c @@ -53,20 +53,25 @@ THIS SOFTWARE. char* #ifdef KR_headers -g_xfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize; +g_xfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize; #else -g_xfmt(char *buf, void *V, int ndig, unsigned bufsize) +g_xfmt(char *buf, void *V, int ndig, size_t bufsize) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 }; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0, Int_max }; char *b, *s, *se; ULong bits[2], sign; UShort *L; int decpt, ex, i, mode; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif if (ndig < 0) ndig = 0; - if (bufsize < (unsigned)(ndig + 10)) + if (bufsize < (size_t)(ndig + 10)) return 0; L = (UShort *)V; @@ -76,14 +81,14 @@ g_xfmt(char *buf, void *V, int ndig, unsigned bufsize) if ( (ex = L[_0] & 0x7fff) !=0) { if (ex == 0x7fff) { /* Infinity or NaN */ - if (bits[0] | bits[1]) - b = strcp(buf, "NaN"); - else { + if (!bits[0] && bits[1]== 0x80000000) { b = buf; if (sign) *b++ = '-'; b = strcp(b, "Infinity"); } + else + b = strcp(buf, "NaN"); return b; } i = STRTOG_Normal; @@ -109,6 +114,6 @@ g_xfmt(char *buf, void *V, int ndig, unsigned bufsize) return 0; mode = 0; } - s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se); - return g__fmt(buf, s, se, decpt, sign); + s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se); + return g__fmt(buf, s, se, decpt, sign, bufsize); } diff --git a/gdtoa/gdtoa.c b/gdtoa/gdtoa.c index 7e8d2b2bf..a4759968a 100644 --- a/gdtoa/gdtoa.c +++ b/gdtoa/gdtoa.c @@ -30,7 +30,6 @@ THIS SOFTWARE. * with " at " changed at "@" and " dot " changed to "."). */ #include "gdtoaimp.h" -#include static Bigint * #ifdef KR_headers @@ -62,7 +61,7 @@ bitstob(ULong *bits, int nbits, int *bbits) *x++ = (*bits >> 16) & ALL_ON; #endif } while(++bits <= be); - i = (int)(x - x0); + i = x - x0; while(!x0[--i]) if (!i) { b->wds = 0; @@ -116,14 +115,15 @@ gdtoa FPI *fpi; int be; ULong *bits; int *kindp, mode, ndigits, *decpt; char **rve; #else - (CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve) + (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve) #endif { /* Arguments ndigits and decpt are similar to the second and third arguments of ecvt and fcvt; trailing zeros are suppressed from the returned string. If not null, *rve is set to point to the end of the return value. If d is +-Infinity or NaN, - then *decpt is set to INT_MAX. + then *decpt is set to 9999. + be = exponent: value = (integer represented by bits) * (2 to the power of be). mode: 0 ==> shortest string that yields d when read in @@ -159,8 +159,8 @@ gdtoa Long L; Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S; double d2, ds; - U d, eps; char *s, *s0; + U d, eps; #ifndef MULTIPLE_THREADS if (dtoa_result) { @@ -177,10 +177,10 @@ gdtoa case STRTOG_Denormal: break; case STRTOG_Infinite: - *decpt = INT_MAX; + *decpt = -32768; return nrv_alloc("Infinity", rve, 8); case STRTOG_NaN: - *decpt = INT_MAX; + *decpt = -32768; return nrv_alloc("NaN", rve, 3); default: return 0; @@ -199,21 +199,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 ); @@ -233,7 +233,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) @@ -248,13 +248,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; } @@ -284,11 +284,14 @@ gdtoa mode -= 4; try_quick = 0; } + else if (i >= -4 - Emin || i < Emin) + 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; @@ -330,10 +333,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; @@ -344,7 +347,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++) @@ -356,30 +359,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; } @@ -388,42 +391,40 @@ 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)) { - while(*--s == '0'){} - s++; - if (dval(d)) + else if (dval(&d) < ds - dval(&eps)) { + if (dval(&d)) inex = STRTOG_Inexlo; - goto ret1; + goto clear_trailing0; } break; } @@ -433,34 +434,34 @@ gdtoa #endif fast_failed: s = s0; - dval(d) = d2; + dval(&d) = d2; k = k0; ilim = ilim0; } /* Do we have a "small" integer? */ - if (be >= 0 && k <= Int_max) { + if (be >= 0 && k <= fpi->int_max) { /* Yes. */ 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) == 0.) + if (dval(&d) == 0.) break; if (i == ilim) { if (rdir) { @@ -469,8 +470,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') @@ -481,8 +487,12 @@ gdtoa } ++*s++; } - else + else { inex = STRTOG_Inexlo; + clear_trailing0: + while(*--s == '0'){} + ++s; + } break; } } @@ -493,13 +503,15 @@ gdtoa m5 = b5; mhi = mlo = 0; if (leftright) { - if (mode < 2) { - i = nbits - bbits; - if (be - i++ < fpi->emin) - /* denormal */ - i = be - fpi->emin + 1; + i = nbits - bbits; + if (be - i++ < fpi->emin && mode != 3 && mode != 5) { + /* denormal */ + i = be - fpi->emin + 1; + if (mode >= 2 && ilim > 0 && ilim < i) + goto small_ilim; } - else { + else if (mode >= 2) { + small_ilim: j = ilim - 1; if (m5 >= j) m5 -= j; @@ -560,28 +572,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 +641,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 +668,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 +717,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') @@ -740,7 +744,7 @@ gdtoa if (b->wds > 1 || b->x[0]) inex = STRTOG_Inexlo; while(*--s == '0'){} - s++; + ++s; } ret: Bfree(S); diff --git a/gdtoa/gdtoa.h b/gdtoa/gdtoa.h index 2707ff7df..49c637beb 100644 --- a/gdtoa/gdtoa.h +++ b/gdtoa/gdtoa.h @@ -66,16 +66,13 @@ THIS SOFTWARE. #else #include "arith.h" #endif +#include /* for size_t */ -/* [RH] On 64-bit Unix systems, long is a 64-bit type. I do not that is - * is what is desired here, so I sacrifice compatibility with systems - * that use 16-bit ints (oh no!) and make Long an int instead. - */ #ifndef Long -typedef int Long; +#define Long int #endif #ifndef ULong -typedef unsigned int ULong; +typedef unsigned Long ULong; #endif #ifndef UShort typedef unsigned short UShort; @@ -111,9 +108,9 @@ typedef unsigned short UShort; /* The following may be or-ed into one of the above values. */ - STRTOG_Neg = 0x08, - STRTOG_Inexlo = 0x10, - STRTOG_Inexhi = 0x20, + STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */ + STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */ + STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */ STRTOG_Inexact = 0x30, STRTOG_Underflow= 0x40, STRTOG_Overflow = 0x80 @@ -126,6 +123,7 @@ FPI { int emax; int rounding; int sudden_underflow; + int int_max; } FPI; enum { /* FPI.rounding values: same as FLT_ROUNDS */ @@ -135,29 +133,31 @@ enum { /* FPI.rounding values: same as FLT_ROUNDS */ FPI_Round_down = 3 }; -// Our strtod is not the CRT's strtod. -#include -#define strtod mystrtod - #ifdef __cplusplus extern "C" { #endif extern char* dtoa ANSI((double d, int mode, int ndigits, int *decpt, int *sign, char **rve)); -extern char* gdtoa ANSI((CONST FPI *fpi, int be, ULong *bits, int *kindp, +extern char* gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)); extern void freedtoa ANSI((char*)); extern float strtof ANSI((CONST char *, char **)); extern double strtod ANSI((CONST char *, char **)); -extern int strtodg ANSI((CONST char*, char**, CONST FPI*, Long*, ULong*)); +extern int strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*)); -extern char* g_ddfmt ANSI((char*, double*, int, unsigned)); -extern char* g_dfmt ANSI((char*, double*, int, unsigned)); -extern char* g_ffmt ANSI((char*, float*, int, unsigned)); -extern char* g_Qfmt ANSI((char*, void*, int, unsigned)); -extern char* g_xfmt ANSI((char*, void*, int, unsigned)); -extern char* g_xLfmt ANSI((char*, void*, int, unsigned)); +extern char* g_ddfmt ANSI((char*, double*, int, size_t)); +extern char* g_ddfmt_p ANSI((char*, double*, int, size_t, int)); +extern char* g_dfmt ANSI((char*, double*, int, size_t)); +extern char* g_dfmt_p ANSI((char*, double*, int, size_t, int)); +extern char* g_ffmt ANSI((char*, float*, int, size_t)); +extern char* g_ffmt_p ANSI((char*, float*, int, size_t, int)); +extern char* g_Qfmt ANSI((char*, void*, int, size_t)); +extern char* g_Qfmt_p ANSI((char*, void*, int, size_t, int)); +extern char* g_xfmt ANSI((char*, void*, int, size_t)); +extern char* g_xfmt_p ANSI((char*, void*, int, size_t, int)); +extern char* g_xLfmt ANSI((char*, void*, int, size_t)); +extern char* g_xLfmt_p ANSI((char*, void*, int, size_t, int)); extern int strtoId ANSI((CONST char*, char**, double*, double*)); extern int strtoIdd ANSI((CONST char*, char**, double*, double*)); diff --git a/gdtoa/gdtoa_fltrnds.h b/gdtoa/gdtoa_fltrnds.h new file mode 100644 index 000000000..33e5f9e53 --- /dev/null +++ b/gdtoa/gdtoa_fltrnds.h @@ -0,0 +1,18 @@ + FPI *fpi, fpi1; + int Rounding; +#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ + Rounding = Flt_Rounds; +#else /*}{*/ + Rounding = 1; + switch(fegetround()) { + case FE_TOWARDZERO: Rounding = 0; break; + case FE_UPWARD: Rounding = 2; break; + case FE_DOWNWARD: Rounding = 3; + } +#endif /*}}*/ + fpi = &fpi0; + if (Rounding != 1) { + fpi1 = fpi0; + fpi = &fpi1; + fpi1.rounding = Rounding; + } diff --git a/gdtoa/gdtoaimp.h b/gdtoa/gdtoaimp.h index f93e03ef5..4dea69e56 100644 --- a/gdtoa/gdtoaimp.h +++ b/gdtoa/gdtoaimp.h @@ -89,12 +89,18 @@ THIS SOFTWARE. * #define IBM for IBM mainframe-style floating-point arithmetic. * #define VAX for VAX-style floating-point arithmetic (D_floating). * #define No_leftright to omit left-right logic in fast floating-point - * computation of dtoa. + * computation of dtoa and gdtoa. This will cause modes 4 and 5 to be + * treated the same as modes 2 and 3 for some inputs. * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3. * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines * that use extended-precision instructions to compute rounded * products and quotients) with IBM. - * #define ROUND_BIASED for IEEE-format with biased rounding. + * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic + * that rounds toward +Infinity. + * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased + * rounding when the underlying floating-point arithmetic uses + * unbiased rounding. This prevent using ordinary floating-point + * arithmetic when the result could be computed with one rounding error. * #define Inaccurate_Divide for IEEE-format with correctly rounded * products but inaccurate quotients, e.g., for Intel i860. * #define NO_LONG_LONG on machines that do not have a "long long" @@ -113,7 +119,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, @@ -126,8 +137,10 @@ THIS SOFTWARE. * conversions of IEEE doubles in single-threaded executions with * 8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with * 4-byte pointers, PRIVATE_MEM >= 7112 appears adequate. - * #define INFNAN_CHECK on IEEE systems to cause strtod to check for - * Infinity and NaN (case insensitively). + * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK + * #defined automatically on IEEE systems. On such systems, + * when INFNAN_CHECK is #defined, strtod checks + * for Infinity and NaN (case insensitively). * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, * strtodg also accepts (case insensitively) strings of the form * NaN(x), where x is a string of hexadecimal digits (optionally @@ -151,7 +164,7 @@ THIS SOFTWARE. * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined. * #define IMPRECISE_INEXACT if you do not care about the setting of * the STRTOG_Inexact bits in the special case of doing IEEE double - * precision conversions (which could also be done by the strtog in + * precision conversions (which could also be done by the strtod in * dtoa.c). * #define NO_HEX_FP to disable recognition of C9x's hexadecimal * floating-point constants. @@ -220,6 +233,10 @@ THIS SOFTWARE. #include "gd_qnan.h" #endif +#ifdef Honor_FLT_ROUNDS +#include +#endif + #ifdef DEBUG #include "stdio.h" #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} @@ -229,13 +246,13 @@ THIS SOFTWARE. #include "string.h" #ifdef KR_headers -#define KR_VOID char +#define Char char #else -#define KR_VOID void +#define Char void #endif #ifdef MALLOC -extern KR_VOID *MALLOC ANSI((size_t)); +extern Char *MALLOC ANSI((size_t)); #else #define MALLOC malloc #endif @@ -307,19 +324,19 @@ extern "C" { #endif #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1 -#error Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. +Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined. #endif typedef union { double d; ULong L[2]; } U; #ifdef IEEE_8087 -#define word0(x) x.L[1] -#define word1(x) x.L[0] +#define word0(x) (x)->L[1] +#define word1(x) (x)->L[0] #else -#define word0(x) x.L[0] -#define word1(x) x.L[1] +#define word0(x) (x)->L[0] +#define word1(x) (x)->L[1] #endif -#define dval(x) x.d +#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 @@ -433,6 +450,11 @@ typedef union { double d; ULong L[2]; } U; #ifndef IEEE_Arith #define ROUND_BIASED +#else +#ifdef ROUND_BIASED_without_Round_Up +#undef ROUND_BIASED +#define ROUND_BIASED +#endif #endif #ifdef RND_PRODQUOT @@ -488,12 +510,12 @@ extern double rnd_prod(double, double), rnd_quot(double, double); #define ALL_ON 0xffff #endif -#ifndef MULTIPLE_THREADS +//#ifndef MULTIPLE_THREADS #define ACQUIRE_DTOA_LOCK(n) /*nothing*/ #define FREE_DTOA_LOCK(n) /*nothing*/ -#endif +//#endif -#define Kmax (sizeof(size_t) << 3) +#define Kmax 9 struct Bigint { @@ -516,12 +538,15 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t)); #define Balloc Balloc_D2A #define Bfree Bfree_D2A +#define InfName InfName_D2A +#define NanName NanName_D2A #define ULtoQ ULtoQ_D2A #define ULtof ULtof_D2A #define ULtod ULtod_D2A #define ULtodd ULtodd_D2A #define ULtox ULtox_D2A #define ULtoxL ULtoxL_D2A +#define add_nanbits add_nanbits_D2A #define any_on any_on_D2A #define b2d b2d_D2A #define bigtens bigtens_D2A @@ -560,9 +585,11 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t)); #define trailz trailz_D2A #define ulp ulp_D2A + extern char *add_nanbits ANSI((char*, size_t, ULong*, int)); extern char *dtoa_result; extern CONST double bigtens[], tens[], tinytens[]; extern unsigned char hexdig[]; + extern const char *InfName[6], *NanName[3]; extern Bigint *Balloc ANSI((int)); extern void Bfree ANSI((Bigint*)); @@ -577,14 +604,14 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t)); extern int cmp ANSI((Bigint*, Bigint*)); extern void copybits ANSI((ULong*, int, Bigint*)); extern Bigint *d2b ANSI((double, int*, int*)); - extern int decrement ANSI((Bigint*)); + extern void decrement ANSI((Bigint*)); extern Bigint *diff ANSI((Bigint*, Bigint*)); extern char *dtoa ANSI((double d, int mode, int ndigits, int *decpt, int *sign, char **rve)); - extern char *g__fmt ANSI((char*, char*, char*, int, ULong)); - extern int gethex ANSI((CONST char**, CONST FPI*, Long*, Bigint**, int)); + extern char *g__fmt ANSI((char*, char*, char*, int, ULong, size_t)); + extern int gethex ANSI((CONST char**, FPI*, Long*, Bigint**, int)); extern void hexdig_init_D2A(Void); - extern int hexnan ANSI((CONST char**, CONST FPI*, ULong*)); + extern int hexnan ANSI((CONST char**, FPI*, ULong*)); extern int hi0bits_D2A ANSI((ULong)); extern Bigint *i2b ANSI((int)); extern Bigint *increment ANSI((Bigint*)); @@ -599,14 +626,14 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t)); extern double ratio ANSI((Bigint*, Bigint*)); extern void rshift ANSI((Bigint*, int)); extern char *rv_alloc ANSI((int)); - extern Bigint *s2b ANSI((CONST char*, int, int, ULong)); + extern Bigint *s2b ANSI((CONST char*, int, int, ULong, int)); extern Bigint *set_ones ANSI((Bigint*, int)); extern char *strcp ANSI((char*, const char*)); - extern int strtoIg ANSI((CONST char*, char**, CONST FPI*, Long*, Bigint**, int*)); + extern int strtoIg ANSI((CONST char*, char**, FPI*, Long*, Bigint**, int*)); 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 } @@ -621,6 +648,10 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t)); * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.) */ #ifdef IEEE_Arith +#ifndef NO_INFNAN_CHECK +#undef INFNAN_CHECK +#define INFNAN_CHECK +#endif #ifdef IEEE_MC68k #define _0 0 #define _1 1 diff --git a/gdtoa/gethex.c b/gdtoa/gethex.c index 7273a6c3b..72da9d326 100644 --- a/gdtoa/gethex.c +++ b/gdtoa/gethex.c @@ -40,22 +40,34 @@ THIS SOFTWARE. gethex(sp, fpi, exp, bp, sign) CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign; #else -gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) +gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign) #endif { Bigint *b; CONST unsigned char *decpt, *s0, *s, *s1; - int esign, havedig, irv, k, n, nbits, up, zret; + int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret; ULong L, lostbits, *x; Long e, e1; #ifdef USE_LOCALE - unsigned char decimalpoint = *localeconv()->decimal_point; + int i; +#ifdef NO_LOCALE_CACHE + const unsigned char *decimalpoint = (unsigned char*)localeconv()->decimal_point; #else -#define decimalpoint '.' + const unsigned char *decimalpoint; + static unsigned char *decimalpoint_cache; + if (!(s0 = decimalpoint_cache)) { + s0 = (unsigned char*)localeconv()->decimal_point; + if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { + strcpy(decimalpoint_cache, s0); + s0 = decimalpoint_cache; + } + } + decimalpoint = s0; +#endif #endif - if (!hexdig['0']) - hexdig_init_D2A(); + /**** if (!hexdig['0']) hexdig_init_D2A(); ****/ + *bp = 0; havedig = 0; s0 = *(CONST unsigned char **)sp + 2; while(s0[havedig] == '0') @@ -65,11 +77,21 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) decpt = 0; zret = 0; e = 0; - if (!hexdig[*s]) { + if (hexdig[*s]) + havedig++; + else { zret = 1; - if (*s != decimalpoint) +#ifdef USE_LOCALE + for(i = 0; decimalpoint[i]; ++i) { + if (s[i] != decimalpoint[i]) + goto pcheck; + } + decpt = s += i; +#else + if (*s != '.') goto pcheck; decpt = ++s; +#endif if (!hexdig[*s]) goto pcheck; while(*s == '0') @@ -81,19 +103,28 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) } while(hexdig[*s]) s++; - if (*s == decimalpoint && !decpt) { +#ifdef USE_LOCALE + if (*s == *decimalpoint && !decpt) { + for(i = 1; decimalpoint[i]; ++i) { + if (s[i] != decimalpoint[i]) + goto pcheck; + } + decpt = s += i; +#else + if (*s == '.' && !decpt) { decpt = ++s; +#endif while(hexdig[*s]) s++; - } + }/*}*/ if (decpt) e = -(((Long)(s-decpt)) << 2); pcheck: s1 = s; + big = esign = 0; switch(*s) { case 'p': case 'P': - esign = 0; switch(*++s) { case '-': esign = 1; @@ -106,29 +137,87 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) break; } e1 = n - 0x10; - while((n = hexdig[*++s]) !=0 && n <= 0x19) + while((n = hexdig[*++s]) !=0 && n <= 0x19) { + if (e1 & 0xf8000000) + big = 1; e1 = 10*e1 + n - 0x10; + } if (esign) e1 = -e1; e += e1; } *sp = (char*)s; - if (zret) { - if (!havedig) - *sp = s0 - 1; + if (!havedig) + *sp = (char*)s0 - 1; + if (zret) return STRTOG_Zero; + if (big) { + if (esign) { + switch(fpi->rounding) { + case FPI_Round_up: + if (sign) + break; + goto ret_tiny; + case FPI_Round_down: + if (!sign) + break; + goto ret_tiny; + } + goto retz; + ret_tiny: + b = Balloc(0); + b->wds = 1; + b->x[0] = 1; + goto dret; + } + switch(fpi->rounding) { + case FPI_Round_near: + goto ovfl1; + case FPI_Round_up: + if (!sign) + goto ovfl1; + goto ret_big; + case FPI_Round_down: + if (sign) + goto ovfl1; + goto ret_big; + } + ret_big: + nbits = fpi->nbits; + n0 = n = nbits >> kshift; + if (nbits & kmask) + ++n; + for(j = n, k = 0; j >>= 1; ++k); + *bp = b = Balloc(k); + b->wds = n; + for(j = 0; j < n0; ++j) + b->x[j] = ALL_ON; + if (n > n0) + b->x[j] = ULbits >> (ULbits - (nbits & kmask)); + *exp = fpi->emin; + return STRTOG_Normal | STRTOG_Inexlo; } - n = (int)(s1 - s0 - 1); - for(k = 0; n > 7; n >>= 1) + n = s1 - s0 - 1; + for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1) k++; b = Balloc(k); x = b->x; n = 0; L = 0; +#ifdef USE_LOCALE + for(i = 0; decimalpoint[i+1]; ++i); +#endif while(s1 > s0) { - if (*--s1 == decimalpoint) +#ifdef USE_LOCALE + if (*--s1 == decimalpoint[i]) { + s1 -= i; continue; - if (n == 32) { + } +#else + if (*--s1 == '.') + continue; +#endif + if (n == ULbits) { *x++ = L; L = 0; n = 0; @@ -137,8 +226,8 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) n += 4; } *x++ = L; - b->wds = n = (int)(x - b->x); - n = 32*n - hi0bits(L); + b->wds = n = x - b->x; + n = ULbits*n - hi0bits(L); nbits = fpi->nbits; lostbits = 0; x = b->x; @@ -149,7 +238,7 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) k = n - 1; if (x[k>>kshift] & 1 << (k & kmask)) { lostbits = 2; - if (k > 1 && any_on(b,k-1)) + if (k > 0 && any_on(b,k)) lostbits = 3; } } @@ -165,7 +254,10 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) if (e > fpi->emax) { ovfl: Bfree(b); - *bp = 0; + ovfl1: +#ifndef NO_ERRNO + errno = ERANGE; +#endif return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; } irv = STRTOG_Normal; @@ -185,15 +277,22 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign) case FPI_Round_down: if (sign) { one_bit: - *exp = fpi->emin; x[0] = b->wds = 1; + dret: *bp = b; + *exp = fpi->emin; +#ifndef NO_ERRNO + errno = ERANGE; +#endif return STRTOG_Denormal | STRTOG_Inexhi | STRTOG_Underflow; } } Bfree(b); - *bp = 0; + retz: +#ifndef NO_ERRNO + errno = ERANGE; +#endif return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow; } k = n - 1; @@ -214,7 +313,7 @@ gethex( CONST char **sp, CONST 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: @@ -233,8 +332,8 @@ gethex( CONST char **sp, CONST 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; diff --git a/gdtoa/gmisc.c b/gdtoa/gmisc.c index 988d3bfbb..8270ef944 100644 --- a/gdtoa/gmisc.c +++ b/gdtoa/gmisc.c @@ -60,7 +60,7 @@ rshift(Bigint *b, int k) while(x < xe) *x1++ = *x++; } - if ((b->wds = (int)(x1 - b->x)) == 0) + if ((b->wds = x1 - b->x) == 0) b->x[0] = 0; } diff --git a/gdtoa/hd_init.c b/gdtoa/hd_init.c index fa6e18dee..d79ae2ec8 100644 --- a/gdtoa/hd_init.c +++ b/gdtoa/hd_init.c @@ -31,6 +31,7 @@ THIS SOFTWARE. #include "gdtoaimp.h" +#if 0 unsigned char hexdig[256]; static void @@ -46,10 +47,31 @@ htinit(unsigned char *h, unsigned char *s, int inc) } void -hexdig_init_D2A(Void) +hexdig_init_D2A(Void) /* Use of hexdig_init omitted 20121220 to avoid a */ + /* race condition when multiple threads are used. */ { #define USC (unsigned char *) htinit(hexdig, USC "0123456789", 0x10); htinit(hexdig, USC "abcdef", 0x10 + 10); htinit(hexdig, USC "ABCDEF", 0x10 + 10); } +#else + unsigned char hexdig[256] = { + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0, + 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, + 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 + }; +#endif diff --git a/gdtoa/hexnan.c b/gdtoa/hexnan.c index 0ee127f56..80721e97a 100644 --- a/gdtoa/hexnan.c +++ b/gdtoa/hexnan.c @@ -54,15 +54,14 @@ L_shift(ULong *x, ULong *x1, int i) hexnan(sp, fpi, x0) CONST char **sp; FPI *fpi; ULong *x0; #else -hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0) +hexnan( CONST char **sp, FPI *fpi, ULong *x0) #endif { ULong c, h, *x, *x1, *xe; CONST char *s; int havedig, hd0, i, nbits; - if (!hexdig['0']) - hexdig_init_D2A(); + /**** if (!hexdig['0']) hexdig_init_D2A(); ****/ nbits = fpi->nbits; x = x0 + (nbits >> kshift); if (nbits & kmask) @@ -72,12 +71,15 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0) havedig = hd0 = i = 0; s = *sp; /* allow optional initial 0x or 0X */ - while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') + while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') { + if (!c) + goto retnan; ++s; + } 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) { @@ -92,8 +94,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0) x1 = x; i = 0; } - while(*(CONST unsigned char*)(s+1) <= ' ') + while((c = *(CONST unsigned char*)(s+1)) <= ' ') { + if (!c) + goto retnan; ++s; + } if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X') && *(CONST unsigned char*)(s+3) > ' ') s += 2; @@ -107,10 +112,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0) do { if (/*(*/ c == ')') { *sp = s + 1; - break; + goto break2; } - } while(c = *++s); + } while((c = *++s)); #endif + retnan: return STRTOG_NaN; } havedig++; @@ -120,8 +126,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0) i = 1; *--x = 0; } - *x = (*x << 4) | h & 0xf; + *x = (*x << 4) | (h & 0xf); } +#ifndef GDTOA_NON_PEDANTIC_NANCHECK + break2: +#endif if (!havedig) return STRTOG_NaN; if (x < x1 && i < 8) diff --git a/gdtoa/misc.c b/gdtoa/misc.c index 51f8f52fe..d13046732 100644 --- a/gdtoa/misc.c +++ b/gdtoa/misc.c @@ -38,11 +38,6 @@ THIS SOFTWARE. #endif #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) static double private_mem[PRIVATE_mem], *pmem_next = private_mem; -#endif - -#ifdef MULTIPLE_THREADS -static void ACQUIRE_DTOA_LOCK(int n); -static void FREE_DTOA_LOCK(int n); #endif Bigint * @@ -60,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 { @@ -70,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; } @@ -94,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); + } } } @@ -109,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) @@ -209,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; @@ -621,8 +626,8 @@ b2d #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; @@ -635,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; @@ -668,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 @@ -679,29 +684,28 @@ 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 int de, k; ULong *x, y, z; - 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 - - dval(d) = _d; + d.d = dd; #ifdef VAX - d0 = word0(d) >> 16 | word0(d) << 16; - d1 = word1(d) >> 16 | word1(d) << 16; + d0 = word0(&d) >> 16 | word0(&d) << 16; + d1 = word1(&d) >> 16 | word1(&d) << 16; #endif #ifdef Pack_32 @@ -725,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 @@ -736,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 @@ -798,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; @@ -851,71 +851,25 @@ 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; } #ifdef NO_STRING_H - KR_VOID * + Char * #ifdef KR_headers -memcpy_D2A(a, b, len) char *a; char *b; size_t len; +memcpy_D2A(a, b, len) Char *a; Char *b; size_t len; #else 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; } #endif /* NO_STRING_H */ - - -#ifdef MULTIPLE_THREADS -#ifdef _WIN32 - -#undef Bias -#define WIN32_LEAN_AND_MEAN -#include - -static CRITICAL_SECTION dtoa_lock[2]; -static int did_init; - -void ACQUIRE_DTOA_LOCK(int n) -{ - if (!did_init) - { - did_init = 1; - InitializeCriticalSection(&dtoa_lock[0]); - InitializeCriticalSection(&dtoa_lock[1]); - } - EnterCriticalSection(&dtoa_lock[n]); -} - -void FREE_DTOA_LOCK(int n) -{ - LeaveCriticalSection(&dtoa_lock[n]); -} - -#else - -#include - -static pthread_mutex_t dtoa_lock[2] = { PTHREAD_MUTEX_INITIALIZER, PTHREAD_MUTEX_INITIALIZER }; - -void ACQUIRE_DTOA_LOCK(int n) -{ - pthread_mutex_lock(&dtoa_lock[n]); -} - -void FREE_DTOA_LOCK(int n) -{ - pthread_mutex_unlock(&dtoa_lock[n]); -} - -#endif -#endif diff --git a/gdtoa/qnan.c b/gdtoa/qnan.c index 118e74921..ea7e8745b 100644 --- a/gdtoa/qnan.c +++ b/gdtoa/qnan.c @@ -55,11 +55,19 @@ typedef unsigned Long Ulong; #ifdef IEEE_8087 #define _0 1 #define _1 0 +#define _3 3 +#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG) +static int perm[4] = { 0, 1, 2, 3 }; +#endif #define HAVE_IEEE #endif #ifdef IEEE_MC68k #define _0 0 #define _1 1 +#define _3 0 +#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG) +static int perm[4] = { 3, 2, 1, 0 }; +#endif #define HAVE_IEEE #endif @@ -79,31 +87,32 @@ main(void) #endif } U; U a, b, c; +#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG) int i; +#endif a.L[0] = b.L[0] = 0x7f800000; c.f = a.f - b.f; - printf("#define f_QNAN 0x%lx\n", UL c.L[0]); + printf("#define f_QNAN 0x%lx\n", UL (c.L[0] & 0x7fffffff)); a.L[_0] = b.L[_0] = 0x7ff00000; a.L[_1] = b.L[_1] = 0; c.d = a.d - b.d; /* quiet NaN */ - printf("#define d_QNAN0 0x%lx\n", UL c.L[0]); - printf("#define d_QNAN1 0x%lx\n", UL c.L[1]); -#ifdef NO_LONG_LONG - for(i = 0; i < 4; i++) - printf("#define ld_QNAN%d 0xffffffff\n", i); - for(i = 0; i < 5; i++) - printf("#define ldus_QNAN%d 0xffff\n", i); -#else - b.D = c.D = a.d; - if (printf("") < 0) - c.D = 37; /* never executed; just defeat optimization */ - a.L[2] = a.L[3] = 0; - a.D = b.D - c.D; - for(i = 0; i < 4; i++) - printf("#define ld_QNAN%d 0x%lx\n", i, UL a.L[i]); - for(i = 0; i < 5; i++) - printf("#define ldus_QNAN%d 0x%x\n", i, a.u[i]); + c.L[_0] &= 0x7fffffff; + printf("#define d_QNAN0 0x%lx\n", UL c.L[_0]); + printf("#define d_QNAN1 0x%lx\n", UL c.L[_1]); +#ifndef NO_LONG_LONG +#ifdef Gen_ld_QNAN + if (sizeof(a.D) >= 16) { + b.D = c.D = a.d; + if (printf("") < 0) + c.D = 37; /* never executed; just defeat optimization */ + a.L[0] = a.L[1] = a.L[2] = a.L[3] = 0; + a.D = b.D - c.D; + a.L[_3] &= 0x7fffffff; + for(i = 0; i < 4; i++) + printf("#define ld_QNAN%d 0x%lx\n", i, UL a.L[perm[i]]); + } +#endif #endif #endif /* HAVE_IEEE */ return 0; diff --git a/gdtoa/smisc.c b/gdtoa/smisc.c index 858024c9e..f4dbafb21 100644 --- a/gdtoa/smisc.c +++ b/gdtoa/smisc.c @@ -34,9 +34,9 @@ THIS SOFTWARE. Bigint * s2b #ifdef KR_headers - (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9; + (s, nd0, nd, y9, dplen) CONST char *s; int dplen, nd0, nd; ULong y9; #else - (CONST char *s, int nd0, int nd, ULong y9) + (CONST char *s, int nd0, int nd, ULong y9, int dplen) #endif { Bigint *b; @@ -60,10 +60,10 @@ s2b s += 9; do b = multadd(b, 10, *s++ - '0'); while(++i < nd0); - s++; + s += dplen; } else - s += 10; + s += dplen + 9; for(; i < nd; i++) b = multadd(b, 10, *s++ - '0'); return b; @@ -80,30 +80,30 @@ ratio 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 diff --git a/gdtoa/strtoIQ.c b/gdtoa/strtoIQ.c index 000fa3645..9ce5120e6 100644 --- a/gdtoa/strtoIQ.c +++ b/gdtoa/strtoIQ.c @@ -38,7 +38,7 @@ strtoIQ(s, sp, a, b) CONST char *s; char **sp; void *a; void *b; strtoIQ(CONST char *s, char **sp, void *a, void *b) #endif { - static CONST FPI fpi = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI }; + static FPI fpi = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI }; Long exp[2]; Bigint *B[2]; int k, rv[2]; diff --git a/gdtoa/strtoId.c b/gdtoa/strtoId.c index d34a278ee..1c97d382d 100644 --- a/gdtoa/strtoId.c +++ b/gdtoa/strtoId.c @@ -38,7 +38,7 @@ strtoId(s, sp, f0, f1) CONST char *s; char **sp; double *f0, *f1; strtoId(CONST char *s, char **sp, double *f0, double *f1) #endif { - static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; Long exp[2]; Bigint *B[2]; int k, rv[2]; diff --git a/gdtoa/strtoIdd.c b/gdtoa/strtoIdd.c index 921597fe7..40b7936bc 100644 --- a/gdtoa/strtoIdd.c +++ b/gdtoa/strtoIdd.c @@ -39,9 +39,9 @@ strtoIdd(CONST char *s, char **sp, double *f0, double *f1) #endif { #ifdef Sudden_Underflow - static CONST FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; + static FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; #else - static CONST FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; + static FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; #endif Long exp[2]; Bigint *B[2]; diff --git a/gdtoa/strtoIf.c b/gdtoa/strtoIf.c index b937ff079..65ecab2e0 100644 --- a/gdtoa/strtoIf.c +++ b/gdtoa/strtoIf.c @@ -38,7 +38,7 @@ strtoIf(s, sp, f0, f1) CONST char *s; char **sp; float *f0, *f1; strtoIf(CONST char *s, char **sp, float *f0, float *f1) #endif { - static CONST FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; + static FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; Long exp[2]; Bigint *B[2]; int k, rv[2]; diff --git a/gdtoa/strtoIg.c b/gdtoa/strtoIg.c index fd181fc9b..6a17760cf 100644 --- a/gdtoa/strtoIg.c +++ b/gdtoa/strtoIg.c @@ -35,7 +35,7 @@ THIS SOFTWARE. #ifdef KR_headers strtoIg(s00, se, fpi, exp, B, rvp) CONST char *s00; char **se; FPI *fpi; Long *exp; Bigint **B; int *rvp; #else -strtoIg(CONST char *s00, char **se, CONST FPI *fpi, Long *exp, Bigint **B, int *rvp) +strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp) #endif { Bigint *b, *b1; @@ -61,16 +61,20 @@ strtoIg(CONST char *s00, char **se, CONST FPI *fpi, Long *exp, Bigint **B, int * if (rv & STRTOG_Inexlo) { swap = 0; b1 = increment(b1); - if (fpi->sudden_underflow - && (rv & STRTOG_Retmask) == STRTOG_Zero) { - b1->x[0] = 0; - b1->x[nw1] = 1L << nb11; - rv1 += STRTOG_Normal - STRTOG_Zero; - rv1 &= ~STRTOG_Underflow; + if ((rv & STRTOG_Retmask) == STRTOG_Zero) { + if (fpi->sudden_underflow) { + b1->x[0] = 0; + b1->x[nw1] = 1L << nb11; + rv1 += STRTOG_Normal - STRTOG_Zero; + rv1 &= ~STRTOG_Underflow; + goto swapcheck; + } + rv1 &= STRTOG_Inexlo | STRTOG_Underflow | STRTOG_Zero; + rv1 |= STRTOG_Inexhi | STRTOG_Denormal; 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); diff --git a/gdtoa/strtoIx.c b/gdtoa/strtoIx.c index 8e972c976..783a631f0 100644 --- a/gdtoa/strtoIx.c +++ b/gdtoa/strtoIx.c @@ -38,7 +38,7 @@ strtoIx(s, sp, a, b) CONST char *s; char **sp; void *a; void *b; strtoIx(CONST char *s, char **sp, void *a, void *b) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; Long exp[2]; Bigint *B[2]; int k, rv[2]; diff --git a/gdtoa/strtoIxL.c b/gdtoa/strtoIxL.c index 86dbc744a..869bfd16f 100644 --- a/gdtoa/strtoIxL.c +++ b/gdtoa/strtoIxL.c @@ -38,7 +38,7 @@ strtoIxL(s, sp, a, b) CONST char *s; char **sp; void *a; void *b; strtoIxL(CONST char *s, char **sp, void *a, void *b) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; Long exp[2]; Bigint *B[2]; int k, rv[2]; diff --git a/gdtoa/strtod.c b/gdtoa/strtod.c index 1f01921a4..3c2230053 100644 --- a/gdtoa/strtod.c +++ b/gdtoa/strtod.c @@ -42,24 +42,43 @@ THIS SOFTWARE. #ifndef NO_IEEE_Scale #define Avoid_Underflow #undef tinytens -/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ +/* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */ /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, - 9007199254740992.e-256 + 9007199254740992.*9007199254740992.e-256 }; #endif #endif -#include - #ifdef Honor_FLT_ROUNDS -#define Rounding rounding #undef Check_FLT_ROUNDS #define Check_FLT_ROUNDS #else #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 @@ -74,21 +93,55 @@ 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, adj; - U rv, rv0, aadj1; + 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 -#ifdef Honor_FLT_ROUNDS - int rounding; -#endif +#ifdef USE_LOCALE /*{{*/ +#ifdef NO_LOCALE_CACHE + char *decimalpoint = localeconv()->decimal_point; + int dplen = strlen(decimalpoint); +#else + char *decimalpoint; + static char *decimalpoint_cache; + static int dplen; + if (!(s0 = decimalpoint_cache)) { + s0 = localeconv()->decimal_point; + if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { + strcpy(decimalpoint_cache, s0); + s0 = decimalpoint_cache; + } + dplen = strlen(s0); + } + decimalpoint = (char*)s0; +#endif /*NO_LOCALE_CACHE*/ +#else /*USE_LOCALE}{*/ +#define dplen 1 +#endif /*USE_LOCALE}}*/ + +#ifdef Honor_FLT_ROUNDS /*{*/ + int Rounding; +#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */ + Rounding = Flt_Rounds; +#else /*}{*/ + Rounding = 1; + switch(fegetround()) { + case FE_TOWARDZERO: Rounding = 0; break; + case FE_UPWARD: Rounding = 2; break; + case FE_DOWNWARD: Rounding = 3; + } +#endif /*}}*/ +#endif /*}*/ - //_control87(_PC_53, _MCW_PC); sign = nz0 = nz = decpt = 0; - dval(rv) = 0.; + dval(&rv) = 0.; for(s = s00;;s++) switch(*s) { case '-': sign = 1; @@ -111,22 +164,18 @@ strtod } break2: if (*s == '0') { -#ifndef NO_HEX_FP +#ifndef NO_HEX_FP /*{*/ { - static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; Long exp; ULong bits[2]; switch(s[1]) { case 'x': case 'X': { -#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) +#ifdef Honor_FLT_ROUNDS FPI fpi1 = fpi; - switch(fegetround()) { - case FE_TOWARDZERO: fpi1.rounding = 0; break; - case FE_UPWARD: fpi1.rounding = 2; break; - case FE_DOWNWARD: fpi1.rounding = 3; - } + fpi1.rounding = Rounding; #else #define fpi1 fpi #endif @@ -146,7 +195,7 @@ strtod goto ret; } } -#endif +#endif /*}*/ nz0 = 1; while(*++s == '0') ; if (!*s) @@ -161,13 +210,17 @@ strtod z = 10*z + c - '0'; nd0 = nd; #ifdef USE_LOCALE - if (c == *localeconv()->decimal_point) + if (c == *decimalpoint) { + for(i = 1; decimalpoint[i]; ++i) + if (s[i] != decimalpoint[i]) + goto dig_done; + s += i; + c = *s; #else - if (c == '.') -#endif - { - decpt = 1; + if (c == '.') { c = *++s; +#endif + decpt = 1; if (!nd) { for(; c == '0'; c = *++s) nz++; @@ -196,7 +249,7 @@ strtod nz = 0; } } - } + }/*}*/ dig_done: e = 0; if (c == 'e' || c == 'E') { @@ -240,7 +293,7 @@ strtod #ifdef INFNAN_CHECK /* Check for Nan and Infinity */ ULong bits[2]; - static CONST FPI fpinan = /* only 52 explicit bits */ + static FPI fpinan = /* only 52 explicit bits */ { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; if (!decpt) switch(c) { @@ -250,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; @@ -262,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 @@ -292,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 @@ -310,6 +363,7 @@ strtod ) { if (!e) goto ret; +#ifndef ROUND_BIASED_without_Round_Up if (e > 0) { if (e <= Ten_pmax) { #ifdef VAX @@ -318,11 +372,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 } @@ -334,25 +388,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; } @@ -362,14 +416,15 @@ 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 +#endif /* ROUND_BIASED_without_Round_Up */ } e1 += nd - k; @@ -383,12 +438,12 @@ strtod scale = 0; #endif #ifdef Honor_FLT_ROUNDS - if ((rounding = Flt_Rounds) >= 2) { + if (Rounding >= 2) { if (sign) - rounding = rounding == 2 ? 0 : 2; + Rounding = Rounding == 2 ? 0 : 2; else - if (rounding != 2) - rounding = 0; + if (Rounding != 2) + Rounding = 0; } #endif #endif /*IEEE_Arith*/ @@ -397,67 +452,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) { + 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; @@ -466,44 +527,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. */ @@ -516,12 +572,12 @@ strtod /* Put digits into bd: true value = bd * 10^e */ - bd0 = s2b(s0, nd0, nd, y); + bd0 = s2b(s0, nd0, nd, y, dplen); 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) { @@ -538,16 +594,23 @@ strtod bd2 -= bbe; bs2 = bb2; #ifdef Honor_FLT_ROUNDS - if (rounding != 1) + if (Rounding != 1) 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 @@ -557,7 +620,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 @@ -596,7 +659,7 @@ strtod delta->sign = 0; i = cmp(delta, bs); #ifdef Honor_FLT_ROUNDS - if (rounding != 1) { + if (Rounding != 1) { if (i < 0) { /* Error is less than an ulp */ if (!delta->x[0] && delta->wds <= 1) { @@ -606,17 +669,17 @@ strtod #endif break; } - if (rounding) { + 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 @@ -625,63 +688,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) { - if (!((rounding>>1) ^ dsign)) + 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)); - if (dsign) - dval(rv) += adj; + dval(&adj) *= ulp(&rv); + if (dsign) { + if (word0(&rv) == Big0 && word1(&rv) == Big1) + goto ovfl; + dval(&rv) += adj.d; + } else - dval(rv) -= adj; + dval(&rv) -= adj.d; goto cont; } #endif /*Honor_FLT_ROUNDS*/ @@ -690,12 +756,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 ) { @@ -720,32 +786,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 @@ -760,7 +828,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 ==> */ @@ -771,10 +839,10 @@ strtod } } #endif /*Avoid_Underflow*/ - L = (word0(rv) & Exp_mask) - Exp_msk1; -#endif /*Sudden_Underflow}*/ - word0(rv) = L | Bndry_mask1; - word1(rv) = 0xffffffff; + L = (word0(&rv) & Exp_mask) - Exp_msk1; +#endif /*Sudden_Underflow}}*/ + word0(&rv) = L | Bndry_mask1; + word1(&rv) = 0xffffffff; #ifdef IBM goto cont; #else @@ -782,16 +850,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 } @@ -803,14 +888,14 @@ strtod } if ((aadj = ratio(delta, bs)) <= 2.) { if (dsign) - aadj = dval(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.; - dval(aadj1) = -1.; + dval(&aadj1) = -1.; } else { /* special case -- power of FLT_RADIX to be */ @@ -820,45 +905,45 @@ strtod aadj = 1./FLT_RADIX; else aadj *= 0.5; - dval(aadj1) = -aadj; + dval(&aadj1) = -aadj; } } else { aadj *= 0.5; - dval(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) - dval(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 = dval(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 @@ -867,58 +952,58 @@ strtod if ((z = (ULong)aadj) <= 0) z = 1; aadj = z; - dval(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 = dval(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) @@ -928,7 +1013,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; } @@ -942,12 +1027,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) @@ -955,32 +1045,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 */ - if (word0(rv) == 0 && word1(rv) == 0) +#ifdef IEEE_Arith + if (!(word0(&rv) & Exp_mask)) +#else + 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); } diff --git a/gdtoa/strtodI.c b/gdtoa/strtodI.c index 992702e25..0b7b8a45c 100644 --- a/gdtoa/strtodI.c +++ b/gdtoa/strtodI.c @@ -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; @@ -55,14 +55,10 @@ strtodI(s, sp, dd) CONST char *s; char **sp; double *dd; strtodI(CONST char *s, char **sp, double *dd) #endif { - static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; 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; } diff --git a/gdtoa/strtodg.c b/gdtoa/strtodg.c index 2d154148d..c2e3365c7 100644 --- a/gdtoa/strtodg.c +++ b/gdtoa/strtodg.c @@ -89,7 +89,7 @@ increment(Bigint *b) return b; } - int + void #ifdef KR_headers decrement(b) Bigint *b; #else @@ -119,7 +119,6 @@ decrement(Bigint *b) *x++ = y & 0xffff; } while(borrow && x < xe); #endif - return STRTOG_Inexlo; } static int @@ -173,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, CONST 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; @@ -183,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) { @@ -206,9 +205,9 @@ rvOK goto ret; } switch(rd) { - case 1: + case 1: /* round down (toward -Infinity) */ goto trunc; - case 2: + case 2: /* round up (toward +Infinity) */ break; default: /* round near */ k = bdif - 1; @@ -292,15 +291,12 @@ rvOK static int #ifdef KR_headers -mantbits(d) double d; +mantbits(d) U *d; #else -mantbits(double _d) +mantbits(U *d) #endif { ULong L; - U d; - - dval(d) = _d; #ifdef VAX L = word1(d) << 16 | word1(d) >> 16; if (L) @@ -322,7 +318,7 @@ strtodg (s00, se, fpi, exp, bits) CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits; #else - (CONST char *s00, char **se, CONST FPI *fpi, Long *exp, ULong *bits) + (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits) #endif { int abe, abits, asub; @@ -332,14 +328,35 @@ strtodg int sudden_underflow; CONST char *s, *s0, *s1; double adj0, tol; - U adj, rv; Long L; - ULong y, z; + U adj, rv; + ULong *b, *be, y, z; Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; +#ifdef USE_LOCALE /*{{*/ +#ifdef NO_LOCALE_CACHE + char *decimalpoint = localeconv()->decimal_point; + int dplen = strlen(decimalpoint); +#else + char *decimalpoint; + static char *decimalpoint_cache; + static int dplen; + if (!(s0 = decimalpoint_cache)) { + s0 = localeconv()->decimal_point; + if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) { + strcpy(decimalpoint_cache, s0); + s0 = decimalpoint_cache; + } + dplen = strlen(s0); + } + decimalpoint = (char*)s0; +#endif /*NO_LOCALE_CACHE*/ +#else /*USE_LOCALE}{*/ +#define dplen 1 +#endif /*USE_LOCALE}}*/ 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) { @@ -394,13 +411,17 @@ strtodg z = 10*z + c - '0'; nd0 = nd; #ifdef USE_LOCALE - if (c == *localeconv()->decimal_point) + if (c == *decimalpoint) { + for(i = 1; decimalpoint[i]; ++i) + if (s[i] != decimalpoint[i]) + goto dig_done; + s += i; + c = *s; #else - if (c == '.') -#endif - { - decpt = 1; + if (c == '.') { c = *++s; +#endif + decpt = 1; if (!nd) { for(; c == '0'; c = *++s) nz++; @@ -429,7 +450,7 @@ strtodg nz = 0; } } - } + }/*}*/ dig_done: e = 0; if (c == 'e' || c == 'E') { @@ -527,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) { @@ -541,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; @@ -556,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; } @@ -595,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 @@ -644,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); @@ -688,7 +709,7 @@ strtodg /* Put digits into bd: true value = bd * 10^e */ - bd0 = s2b(s0, nd0, nd, y); + bd0 = s2b(s0, nd0, nd, y, dplen); for(;;) { bd = Balloc(bd0->k); @@ -822,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); @@ -839,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) { @@ -853,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 = (Long)adj0; adj0 -= L; switch(rd) { @@ -880,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) { @@ -893,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) @@ -947,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; } @@ -980,6 +1001,29 @@ strtodg Bfree(bd0); Bfree(delta); if (rve > fpi->emax) { + switch(fpi->rounding & 3) { + case FPI_Round_near: + goto huge; + case FPI_Round_up: + if (!sign) + goto huge; + break; + case FPI_Round_down: + if (sign) + goto huge; + } + /* Round to largest representable magnitude */ + Bfree(rvb); + rvb = 0; + irv = STRTOG_Normal | STRTOG_Inexlo; + *exp = fpi->emax; + b = bits; + be = b + ((fpi->nbits + 31) >> 5); + while(b < be) + *b++ = -1; + if ((j = fpi->nbits & 0x1f)) + *--be >>= (32 - j); + goto ret; huge: rvb->wds = 0; irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; @@ -994,12 +1038,19 @@ strtodg if (sudden_underflow) { rvb->wds = 0; irv = STRTOG_Underflow | STRTOG_Inexlo; +#ifndef NO_ERRNO + errno = ERANGE; +#endif } else { irv = (irv & ~STRTOG_Retmask) | (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); - if (irv & STRTOG_Inexact) + if (irv & STRTOG_Inexact) { irv |= STRTOG_Underflow; +#ifndef NO_ERRNO + errno = ERANGE; +#endif + } } } if (se) diff --git a/gdtoa/strtodnrp.c b/gdtoa/strtodnrp.c index 3d0297f52..19a769f0b 100644 --- a/gdtoa/strtodnrp.c +++ b/gdtoa/strtodnrp.c @@ -44,7 +44,7 @@ strtod(s, sp) CONST char *s; char **sp; strtod(CONST char *s, char **sp) #endif { - static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; ULong bits[2]; Long exp; int k; diff --git a/gdtoa/strtof.c b/gdtoa/strtof.c index 561264a1d..a8beb3520 100644 --- a/gdtoa/strtof.c +++ b/gdtoa/strtof.c @@ -38,13 +38,18 @@ strtof(s, sp) CONST char *s; char **sp; strtof(CONST char *s, char **sp) #endif { - static CONST FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; + static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; ULong bits[1]; Long exp; int k; union { ULong L[1]; float f; } u; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); switch(k & STRTOG_Retmask) { case STRTOG_NoNumber: case STRTOG_Zero: @@ -53,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: diff --git a/gdtoa/strtopQ.c b/gdtoa/strtopQ.c index 9eefd768f..2acf7e910 100644 --- a/gdtoa/strtopQ.c +++ b/gdtoa/strtopQ.c @@ -49,6 +49,9 @@ THIS SOFTWARE. #define _3 0 #endif + extern ULong NanDflt_Q_D2A[4]; + + int #ifdef KR_headers strtopQ(s, sp, V) CONST char *s; char **sp; void *V; @@ -56,13 +59,18 @@ strtopQ(s, sp, V) CONST char *s; char **sp; void *V; strtopQ(CONST char *s, char **sp, void *V) #endif { - static CONST FPI fpi = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI }; + static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI }; ULong bits[4]; Long exp; int k; ULong *L = (ULong*)V; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); switch(k & STRTOG_Retmask) { case STRTOG_NoNumber: case STRTOG_Zero: @@ -90,10 +98,10 @@ strtopQ(CONST char *s, char **sp, void *V) break; case STRTOG_NaN: - L[0] = ld_QNAN0; - L[1] = ld_QNAN1; - L[2] = ld_QNAN2; - L[3] = ld_QNAN3; + L[_0] = NanDflt_Q_D2A[3]; + L[_1] = NanDflt_Q_D2A[2]; + L[_2] = NanDflt_Q_D2A[1]; + L[_3] = NanDflt_Q_D2A[0]; } if (k & STRTOG_Neg) L[_0] |= 0x80000000L; diff --git a/gdtoa/strtopd.c b/gdtoa/strtopd.c index 4c8bbaf42..0fb35daea 100644 --- a/gdtoa/strtopd.c +++ b/gdtoa/strtopd.c @@ -38,12 +38,17 @@ strtopd(s, sp, d) char *s; char **sp; double *d; strtopd(CONST char *s, char **sp, double *d) #endif { - static CONST FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; ULong bits[2]; Long exp; int k; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi0, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); ULtod((ULong*)d, bits, exp, k); return k; } diff --git a/gdtoa/strtopdd.c b/gdtoa/strtopdd.c index d2c6d9b89..738372d88 100644 --- a/gdtoa/strtopdd.c +++ b/gdtoa/strtopdd.c @@ -39,9 +39,9 @@ strtopdd(CONST char *s, char **sp, double *dd) #endif { #ifdef Sudden_Underflow - static CONST FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; + static FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; #else - static CONST FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; + static FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; #endif ULong bits[4]; Long exp; @@ -51,8 +51,13 @@ strtopdd(CONST char *s, char **sp, double *dd) ULong L[4]; } U; U *u; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - rv = strtodg(s, sp, &fpi, &exp, bits); + rv = strtodg(s, sp, fpi, &exp, bits); u = (U*)dd; switch(rv & STRTOG_Retmask) { case STRTOG_NoNumber: @@ -62,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; @@ -74,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; } } @@ -87,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; } } @@ -100,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: @@ -119,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; @@ -131,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: diff --git a/gdtoa/strtopf.c b/gdtoa/strtopf.c index 806e14f6f..23ca5cbe5 100644 --- a/gdtoa/strtopf.c +++ b/gdtoa/strtopf.c @@ -38,12 +38,17 @@ strtopf(s, sp, f) CONST char *s; char **sp; float *f; strtopf(CONST char *s, char **sp, float *f) #endif { - static CONST FPI fpi = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; + static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; ULong bits[1], *L; Long exp; int k; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); L = (ULong*)f; switch(k & STRTOG_Retmask) { case STRTOG_NoNumber: @@ -53,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: diff --git a/gdtoa/strtopx.c b/gdtoa/strtopx.c index ecd0b3991..32192c572 100644 --- a/gdtoa/strtopx.c +++ b/gdtoa/strtopx.c @@ -31,6 +31,8 @@ THIS SOFTWARE. #include "gdtoaimp.h" + extern UShort NanDflt_ldus_D2A[5]; + #undef _0 #undef _1 @@ -58,13 +60,18 @@ strtopx(s, sp, V) CONST char *s; char **sp; void *V; strtopx(CONST char *s, char **sp, void *V) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; ULong bits[2]; Long exp; int k; UShort *L = (UShort*)V; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); switch(k & STRTOG_Retmask) { case STRTOG_NoNumber: case STRTOG_Zero: @@ -87,15 +94,16 @@ 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: - L[0] = ldus_QNAN0; - L[1] = ldus_QNAN1; - L[2] = ldus_QNAN2; - L[3] = ldus_QNAN3; - L[4] = ldus_QNAN4; + L[_4] = NanDflt_ldus_D2A[0]; + L[_3] = NanDflt_ldus_D2A[1]; + L[_2] = NanDflt_ldus_D2A[2]; + L[_1] = NanDflt_ldus_D2A[3]; + L[_0] = NanDflt_ldus_D2A[4]; } if (k & STRTOG_Neg) L[_0] |= 0x8000; diff --git a/gdtoa/strtopxL.c b/gdtoa/strtopxL.c index 00810fa35..6166c1e62 100644 --- a/gdtoa/strtopxL.c +++ b/gdtoa/strtopxL.c @@ -31,6 +31,8 @@ THIS SOFTWARE. #include "gdtoaimp.h" + extern ULong NanDflt_xL_D2A[3]; + #undef _0 #undef _1 @@ -54,13 +56,18 @@ strtopxL(s, sp, V) CONST char *s; char **sp; void *V; strtopxL(CONST char *s, char **sp, void *V) #endif { - static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; ULong bits[2]; Long exp; int k; ULong *L = (ULong*)V; +#ifdef Honor_FLT_ROUNDS +#include "gdtoa_fltrnds.h" +#else +#define fpi &fpi0 +#endif - k = strtodg(s, sp, &fpi, &exp, bits); + k = strtodg(s, sp, fpi, &exp, bits); switch(k & STRTOG_Retmask) { case STRTOG_NoNumber: case STRTOG_Zero: @@ -77,13 +84,14 @@ 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: - L[0] = ld_QNAN0; - L[1] = ld_QNAN1; - L[2] = ld_QNAN2; + L[_0] = NanDflt_xL_D2A[2]; + L[_1] = NanDflt_xL_D2A[1]; + L[_2] = NanDflt_xL_D2A[0]; } if (k & STRTOG_Neg) L[_0] |= 0x80000000L; diff --git a/gdtoa/strtorQ.c b/gdtoa/strtorQ.c index 1cb157890..f5fd7bba9 100644 --- a/gdtoa/strtorQ.c +++ b/gdtoa/strtorQ.c @@ -49,6 +49,8 @@ THIS SOFTWARE. #define _3 0 #endif + extern ULong NanDflt_Q_D2A[4]; + void #ifdef KR_headers ULtoQ(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k; @@ -83,10 +85,10 @@ ULtoQ(ULong *L, ULong *bits, Long exp, int k) break; case STRTOG_NaN: - L[0] = ld_QNAN0; - L[1] = ld_QNAN1; - L[2] = ld_QNAN2; - L[3] = ld_QNAN3; + L[_0] = NanDflt_Q_D2A[3]; + L[_1] = NanDflt_Q_D2A[2]; + L[_2] = NanDflt_Q_D2A[1]; + L[_3] = NanDflt_Q_D2A[0]; } if (k & STRTOG_Neg) L[_0] |= 0x80000000L; @@ -99,9 +101,8 @@ strtorQ(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L; strtorQ(CONST char *s, char **sp, int rounding, void *L) #endif { - static CONST FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI }; - CONST FPI *fpi; - FPI fpi1; + static FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI }; + FPI *fpi, fpi1; ULong bits[4]; Long exp; int k; diff --git a/gdtoa/strtord.c b/gdtoa/strtord.c index 1aa7aefe3..dd0769698 100644 --- a/gdtoa/strtord.c +++ b/gdtoa/strtord.c @@ -31,6 +31,8 @@ THIS SOFTWARE. #include "gdtoaimp.h" + extern ULong NanDflt_d_D2A[2]; + void #ifdef KR_headers ULtod(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k; @@ -61,8 +63,8 @@ ULtod(ULong *L, ULong *bits, Long exp, int k) break; case STRTOG_NaN: - L[0] = d_QNAN0; - L[1] = d_QNAN1; + L[_0] = NanDflt_d_D2A[1]; + L[_1] = NanDflt_d_D2A[0]; } if (k & STRTOG_Neg) L[_0] |= 0x80000000L; @@ -75,9 +77,8 @@ strtord(s, sp, rounding, d) CONST char *s; char **sp; int rounding; double *d; strtord(CONST char *s, char **sp, int rounding, double *d) #endif { - static CONST FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; - CONST FPI *fpi; - FPI fpi1; + static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; + FPI *fpi, fpi1; ULong bits[2]; Long exp; int k; diff --git a/gdtoa/strtordd.c b/gdtoa/strtordd.c index 5b505f1b8..62152dbd4 100644 --- a/gdtoa/strtordd.c +++ b/gdtoa/strtordd.c @@ -31,6 +31,8 @@ THIS SOFTWARE. #include "gdtoaimp.h" + extern ULong NanDflt_d_D2A[2]; + void #ifdef KR_headers ULtodd(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k; @@ -48,8 +50,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 +62,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 +75,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 +88,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 +107,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 +119,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: @@ -153,16 +155,17 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k) break; case STRTOG_NaN: - L[0] = L[2] = d_QNAN0; - L[1] = L[3] = d_QNAN1; + L[_0] = L[_0+2] = NanDflt_d_D2A[1]; + L[_1] = L[_1+2] = NanDflt_d_D2A[0]; break; case STRTOG_NaNbits: - L[_1] = (bits[1] >> 21 | bits[2] << 11) & (ULong)0xffffffffL; - L[_0] = bits[2] >> 21 | bits[3] << 11 - | (ULong)0x7ff00000L; - L[2+_1] = bits[0]; - L[2+_0] = bits[1] | (ULong)0x7ff00000L; + L[_1] = (bits[1] >> 20 | bits[2] << 12) & (ULong)0xffffffffL; + L[_0] = bits[2] >> 20 | bits[3] << 12; + L[_0] |= (L[_1] | L[_0]) ? (ULong)0x7ff00000L : (ULong)0x7ff80000L; + L[2+_1] = bits[0] & (ULong)0xffffffffL; + L[2+_0] = bits[1] & 0xfffffL; + L[2+_0] |= (L[2+_1] | L[2+_0]) ? (ULong)0x7ff00000L : (ULong)0x7ff80000L; } if (k & STRTOG_Neg) { L[_0] |= 0x80000000L; @@ -178,12 +181,11 @@ strtordd(CONST char *s, char **sp, int rounding, double *dd) #endif { #ifdef Sudden_Underflow - static CONST FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; + static FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 }; #else - static CONST FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; + static FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 }; #endif - CONST FPI *fpi; - FPI fpi1; + FPI *fpi, fpi1; ULong bits[4]; Long exp; int k; diff --git a/gdtoa/strtorf.c b/gdtoa/strtorf.c index 3892d68d6..99b4ab710 100644 --- a/gdtoa/strtorf.c +++ b/gdtoa/strtorf.c @@ -31,6 +31,8 @@ THIS SOFTWARE. #include "gdtoaimp.h" + extern ULong NanDflt_f_D2A[1]; + void #ifdef KR_headers ULtof(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k; @@ -46,7 +48,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: @@ -58,7 +60,7 @@ ULtof(ULong *L, ULong *bits, Long exp, int k) break; case STRTOG_NaN: - L[0] = f_QNAN; + L[0] = NanDflt_f_D2A[0]; } if (k & STRTOG_Neg) L[0] |= 0x80000000L; @@ -71,9 +73,8 @@ strtorf(s, sp, rounding, f) CONST char *s; char **sp; int rounding; float *f; strtorf(CONST char *s, char **sp, int rounding, float *f) #endif { - static CONST FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; - CONST FPI *fpi; - FPI fpi1; + static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI }; + FPI *fpi, fpi1; ULong bits[1]; Long exp; int k; diff --git a/gdtoa/strtorx.c b/gdtoa/strtorx.c index cc3c0700f..994ce8e63 100644 --- a/gdtoa/strtorx.c +++ b/gdtoa/strtorx.c @@ -51,6 +51,8 @@ THIS SOFTWARE. #define _4 0 #endif + extern UShort NanDflt_ldus_D2A[5]; + void #ifdef KR_headers ULtox(L, bits, exp, k) UShort *L; ULong *bits; Long exp; int k; @@ -80,15 +82,16 @@ 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: - L[0] = ldus_QNAN0; - L[1] = ldus_QNAN1; - L[2] = ldus_QNAN2; - L[3] = ldus_QNAN3; - L[4] = ldus_QNAN4; + L[_4] = NanDflt_ldus_D2A[0]; + L[_3] = NanDflt_ldus_D2A[1]; + L[_2] = NanDflt_ldus_D2A[2]; + L[_1] = NanDflt_ldus_D2A[3]; + L[_0] = NanDflt_ldus_D2A[4]; } if (k & STRTOG_Neg) L[_0] |= 0x8000; @@ -101,9 +104,8 @@ strtorx(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L; strtorx(CONST char *s, char **sp, int rounding, void *L) #endif { - static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; - CONST FPI *fpi; - FPI fpi1; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + FPI *fpi, fpi1; ULong bits[2]; Long exp; int k; diff --git a/gdtoa/strtorxL.c b/gdtoa/strtorxL.c index ffcafa907..bac4a0bb1 100644 --- a/gdtoa/strtorxL.c +++ b/gdtoa/strtorxL.c @@ -47,6 +47,8 @@ THIS SOFTWARE. #define _2 0 #endif + extern ULong NanDflt_xL_D2A[3]; + void #ifdef KR_headers ULtoxL(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k; @@ -69,14 +71,15 @@ ULtoxL(ULong *L, ULong *bits, Long exp, int k) break; case STRTOG_Infinite: - L[_0] = 0x7fff << 16; - L[_1] = L[_2] = 0; + L[_0] = 0x7fff0000; + L[_1] = 0x80000000; + L[_2] = 0; break; case STRTOG_NaN: - L[0] = ld_QNAN0; - L[1] = ld_QNAN1; - L[2] = ld_QNAN2; + L[_0] = NanDflt_xL_D2A[2]; + L[_1] = NanDflt_xL_D2A[1]; + L[_2] = NanDflt_xL_D2A[0]; } if (k & STRTOG_Neg) L[_0] |= 0x80000000L; @@ -89,9 +92,8 @@ strtorxL(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L; strtorxL(CONST char *s, char **sp, int rounding, void *L) #endif { - static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; - CONST FPI *fpi; - FPI fpi1; + static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI }; + FPI *fpi, fpi1; ULong bits[2]; Long exp; int k; diff --git a/gdtoa/ulp.c b/gdtoa/ulp.c index 27639d19e..17e9f862c 100644 --- a/gdtoa/ulp.c +++ b/gdtoa/ulp.c @@ -34,15 +34,14 @@ THIS SOFTWARE. double ulp #ifdef KR_headers - (x) double x; + (x) U *x; #else - (double _x) + (U *x) #endif { Long L; - U x, a; + U a; - dval(x) = _x; L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; #ifndef Sudden_Underflow if (L > 0) { @@ -50,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 dval(a); + return dval(&a); }