2009-04-14 Ozkan Sezer <sezeroz@gmail.com>

* gdtoa/dmisc.c, gdtoa/dtoa.c, gdtoa/gdtoa.c, gdtoa/gdtoa.h,
	gdtoa/gdtoaimp.h, gdtoa/misc.c, gdtoa/smisc.c, gdtoa/strtodg.c,
	gdtoa/strtof.c, gdtoa/ulp.c: merged the aliasing violation changes from
	the latest netlib.org sources.


git-svn-id: svn+ssh://svn.code.sf.net/p/mingw-w64/code/trunk@769 4407c894-4637-0410-b4f5-ada5f102cad1
This commit is contained in:
Ozkan Sezer 2009-04-14 11:11:44 +00:00
parent 987c9aad6a
commit c79f961db6
11 changed files with 223 additions and 226 deletions

View File

@ -1,3 +1,10 @@
2009-04-14 Ozkan Sezer <sezeroz@gmail.com>
* gdtoa/dmisc.c, gdtoa/dtoa.c, gdtoa/gdtoa.c, gdtoa/gdtoa.h,
gdtoa/gdtoaimp.h, gdtoa/misc.c, gdtoa/smisc.c, gdtoa/strtodg.c,
gdtoa/strtof.c, gdtoa/ulp.c: merged the aliasing violation changes from
the latest netlib.org sources.
2009-04-13 Ozkan Sezer <sezeroz@gmail.com>
* gdtoa/dtoa.c, gdtoa/gdtoa.c, gdtoa/gdtoaimp.h, gdtoa/misc.c,

View File

@ -31,6 +31,10 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
#ifndef MULTIPLE_THREADS
char *dtoa_result;
#endif
char *rv_alloc (int i)
{
int j, k, *r;

View File

@ -109,7 +109,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
spec_case, try_quick;
Long L;
@ -118,17 +118,16 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
ULong x;
#endif
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
union _dbl_union d, d2, eps;
double ds;
char *s, *s0;
#ifdef Honor_FLT_ROUNDS
int rounding;
#endif
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
union _dbl_union d, d2, eps;
#ifdef Honor_FLT_ROUNDS
int rounding;
#endif
d.d = val;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
__freedtoa(dtoa_result);
@ -136,23 +135,24 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
}
#endif
if (word0(d) & Sign_bit) {
d.d = val;
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 ((word0(d) & Exp_mask) == Exp_mask)
if ((word0(&d) & Exp_mask) == Exp_mask)
{
/* Infinity or NaN */
*decpt = 9999;
if (!word1(d) && !(word0(d) & 0xfffff))
if (!word1(&d) && !(word0(&d) & 0xfffff))
return nrv_alloc("Infinity", rve, 8);
return nrv_alloc("NaN", rve, 3);
}
if (!dval(d)) {
if (!dval(&d)) {
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
@ -171,22 +171,22 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
}
#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;
/* 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 );
@ -211,21 +211,21 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
/* 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;
}
@ -264,10 +264,11 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
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;
@ -299,7 +300,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
if (ilim >= 0 && ilim <= Quick_max && try_quick) {
/* 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 */
@ -309,7 +310,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
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++)
@ -317,32 +318,32 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
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;
}
@ -351,34 +352,34 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = 0.5/tens[ilim-1] - dval(eps);
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = dval(d);
dval(d) -= L;
L = dval(&d);
dval(&d) -= L;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps))
if (dval(&d) < dval(&eps))
goto ret1;
if (1. - dval(d) < dval(eps))
if (1. - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d));
if (!(dval(d) -= L))
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d));
if (!(dval(&d) -= L))
ilim = i;
*s++ = '0' + (int)L;
if (i == ilim) {
if (dval(d) > 0.5 + dval(eps))
if (dval(&d) > 0.5 + dval(&eps))
goto bump_up;
else if (dval(d) < 0.5 - dval(eps)) {
else if (dval(&d) < 0.5 - dval(&eps)) {
while(*--s == '0');
s++;
goto ret1;
@ -391,7 +392,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
#endif
fast_failed:
s = s0;
dval(d) = dval(d2);
dval(&d) = dval(&d2);
k = k0;
ilim = ilim0;
}
@ -403,22 +404,22 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
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
@ -432,8 +433,8 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
case 2: goto bump_up;
}
#endif
dval(d) += dval(d);
if (dval(d) > ds || (dval(d) == ds && L & 1)) {
dval(&d) += dval(&d);
if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {
bump_up:
while(*--s == '9')
if (s == s0) {
@ -494,9 +495,9 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
&& 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 */
@ -582,7 +583,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(d) & 1)
if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
#ifdef Honor_FLT_ROUNDS
&& rounding >= 1
#endif
@ -601,7 +602,7 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
#endif
if (j < 0 || (j == 0 && mode != 1
#ifndef ROUND_BIASED
&& !(word1(d) & 1)
&& !(word1(&d) & 1)
#endif
)) {
if (!b->x[0] && b->wds <= 1) {
@ -706,9 +707,9 @@ char *__dtoa (double val, int mode, int ndigits, int *decpt, int *sign, char **r
#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)

View File

@ -140,14 +140,14 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be0, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0, inex;
int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
int rdir, s2, s5, spec_case, try_quick;
Long L;
Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
double d2, ds;
union _dbl_union d, eps;
char *s, *s0;
union _dbl_union d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
@ -186,17 +186,17 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
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;
/* 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 );
@ -212,7 +212,7 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
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)
@ -224,9 +224,9 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
word0(d) += (be + bbits - 1) << Exp_shift;
word0(&d) += (be + bbits - 1) << Exp_shift;
if (k >= 0 && k <= Ten_pmax) {
if (dval(d) < tens[k])
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
@ -257,10 +257,11 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
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;
@ -302,7 +303,7 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
/* Try to get by with floating-point arithmetic. */
i = 0;
d2 = dval(d);
d2 = dval(&d);
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
@ -312,7 +313,7 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
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++)
@ -324,30 +325,30 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
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;
}
@ -356,40 +357,40 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
/* 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)) {
else if (dval(&d) < ds - dval(&eps)) {
while(*--s == '0'){}
s++;
if (dval(d))
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
@ -401,7 +402,7 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
#endif
fast_failed:
s = s0;
dval(d) = d2;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
@ -413,22 +414,22 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = dval(d) / ds;
dval(d) -= L*ds;
for(i = 1;; i++, dval(&d) *= 10.) {
L = dval(&d) / ds;
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (dval(d) == 0.)
if (dval(&d) == 0.)
break;
if (i == ilim) {
if (rdir) {
@ -437,8 +438,8 @@ char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
inex = STRTOG_Inexlo;
goto ret1;
}
dval(d) += dval(d);
if (dval(d) > ds || (dval(d) == ds && L & 1)) {
dval(&d) += dval(&d);
if (dval(&d) > ds || (dval(&d) == ds && L & 1)) {
bump_up:
inex = STRTOG_Inexhi;
while(*--s == '9')

View File

@ -98,10 +98,10 @@ extern char* __gdtoa (FPI *fpi, int be, ULong *bits, int *kindp,
int mode, int ndigits, int *decpt, char **rve);
extern void __freedtoa (char *);
extern int __strtodg (const char *, char **, FPI *, Long *, ULong *);
extern float __strtof (const char *, char **);
extern double __strtod (const char *, char **);
extern long double strtold (const char *, char **);
extern int __strtodg (const char *, char **, FPI *, Long *, ULong *);
extern char* __g__fmt (char*, char*, char*, int, ULong);
extern char* __g_dfmt (char*, double*, int, unsigned);

View File

@ -157,11 +157,6 @@ THIS SOFTWARE.
* #define NO_STRING_H to use private versions of memcpy.
* On some K&R systems, it may also be necessary to
* #define DECLARE_SIZE_T in this case.
* #define YES_ALIAS to permit aliasing certain double values with
* arrays of ULongs. This leads to slightly better code with
* some compilers and was always used prior to 19990916, but it
* is not strictly legal and can cause trouble with aggressively
* optimizing compilers (e.g., gcc 2.95.1 under -O2).
* #define USE_LOCALE to use the current locale's decimal_point value.
*/
@ -261,18 +256,14 @@ Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
typedef union _dbl_union { double d; ULong L[2]; } dbl_union;
/* gcc >= 4.4 at -O2 and higher is seriously allergic to aliasing
violations, unions must be used instead of typecast assignment:
in the following macros, the x argument must be of 'dbl_union'
type. */
#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
@ -567,7 +558,7 @@ extern Bigint *set_ones (Bigint*, int);
extern char *strcp (char*, const char*);
extern Bigint *sum (Bigint*, Bigint*);
extern int trailz (Bigint*);
extern double ulp (double);
extern double ulp (dbl_union *);
#ifdef __cplusplus
}

View File

@ -30,16 +30,16 @@ THIS SOFTWARE.
* with " at " changed at "@" and " dot " changed to "."). */
#ifdef __MINGW32__
/* we have to include windows.h before gdtoa headers, otherwise
defines cause conflicts. */
#if defined(__MINGW32__) || defined(__MINGW64__)
/* we have to include windows.h before gdtoa
headers, otherwise defines cause conflicts. */
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#define NLOCKS 2
#ifdef USE_WIN32_SL
/* Use spin locks. */
/* Use spin locks. */
static long dtoa_sl[NLOCKS];
#define ACQUIRE_DTOA_LOCK(n) \
@ -47,7 +47,8 @@ static long dtoa_sl[NLOCKS];
Sleep (0);
#define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0);
#else
#else /* USE_WIN32_SL */
#include <stdlib.h>
static CRITICAL_SECTION dtoa_CritSec[NLOCKS];
static long dtoa_CS_init = 0;
@ -101,16 +102,12 @@ static void dtoa_unlock (int n)
#define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n)
#define FREE_DTOA_LOCK(n) dtoa_unlock(n)
#endif
#endif /* USE_WIN32_SL */
#endif /* __MINGW32__ */
#endif /* __MINGW32__ / __MINGW64__ */
#include "gdtoaimp.h"
#ifndef MULTIPLE_THREADS
char *dtoa_result;
#endif
static Bigint *freelist[Kmax+1];
#ifndef Omit_Private_Memory
#ifndef PRIVATE_MEM
@ -630,52 +627,53 @@ double b2d (Bigint *a, int *e)
*e = 32 - k;
#ifdef Pack_32
if (k < Ebits) {
word0(d) = Exp_1 | y >> (Ebits - k);
word0(&d) = Exp_1 | y >> (Ebits - k);
w = xa > xa0 ? *--xa : 0;
word1(d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
if (k -= Ebits) {
word0(d) = Exp_1 | y << k | z >> (32 - k);
word0(&d) = Exp_1 | y << k | z >> (32 - k);
y = xa > xa0 ? *--xa : 0;
word1(d) = z << k | y >> (32 - k);
word1(&d) = z << k | y >> (32 - k);
}
else {
word0(d) = Exp_1 | y;
word1(d) = z;
word0(&d) = Exp_1 | y;
word1(&d) = z;
}
#else
if (k < Ebits + 16) {
z = xa > xa0 ? *--xa : 0;
word0(d) = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
word0(&d) = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
w = xa > xa0 ? *--xa : 0;
y = xa > xa0 ? *--xa : 0;
word1(d) = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
word1(&d) = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
w = xa > xa0 ? *--xa : 0;
k -= Ebits + 16;
word0(d) = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
word0(&d) = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
y = xa > xa0 ? *--xa : 0;
word1(d) = w << k + 16 | y << k;
word1(&d) = w << k + 16 | y << k;
#endif
ret_d:
return dval(d);
return dval(&d);
}
Bigint *d2b (double val, int *e, int *bits)
{
Bigint *b;
union _dbl_union d;
#ifndef Sudden_Underflow
int i;
#endif
int de, k;
ULong *x, y, z;
union _dbl_union d;
d.d = val;
#ifdef Pack_32
b = Balloc(1);
#else
@ -683,17 +681,17 @@ Bigint *d2b (double val, int *e, int *bits)
#endif
x = b->x;
z = word0(d) & Frac_mask;
word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
z = word0(&d) & Frac_mask;
word0(&d) &= 0x7fffffff; /* clear sign bit, which we ignore */
#ifdef Sudden_Underflow
de = (int)(word0(d) >> Exp_shift);
de = (int)(word0(&d) >> Exp_shift);
z |= Exp_msk11;
#else
if ( (de = (int)(word0(d) >> Exp_shift)) !=0)
if ( (de = (int)(word0(&d) >> Exp_shift)) !=0)
z |= Exp_msk1;
#endif
#ifdef Pack_32
if ( (y = word1(d)) !=0) {
if ( (y = word1(&d)) !=0) {
if ( (k = lo0bits(&y)) !=0) {
x[0] = y | z << (32 - k);
z >>= k;
@ -719,7 +717,7 @@ Bigint *d2b (double val, int *e, int *bits)
k += 32;
}
#else
if ( (y = word1(d)) !=0) {
if ( (y = word1(&d)) !=0) {
if ( (k = lo0bits(&y)) !=0)
if (k >= 16) {
x[0] = y | z << 32 - k & 0xffff;

View File

@ -68,16 +68,16 @@ double ratio (Bigint *a, Bigint *b)
union _dbl_union 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);
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;
}
return dval(da) / dval(db);
return dval(&da) / dval(&db);
}
#ifdef INFNAN_CHECK

View File

@ -145,7 +145,7 @@ Bigint *set_ones (Bigint *b, int n)
return b;
}
static int rvOK (double d, FPI *fpi, Long *exp, ULong *bits,
static int rvOK (dbl_union *d, FPI *fpi, Long *exp, ULong *bits,
int exact, int rd, int *irv)
{
Bigint *b;
@ -153,7 +153,7 @@ static int rvOK (double d, FPI *fpi, Long *exp, ULong *bits,
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) {
@ -260,12 +260,9 @@ static int rvOK (double d, FPI *fpi, Long *exp, ULong *bits,
return rv;
}
static int mantbits (double val)
static int mantbits (dbl_union *d)
{
ULong L;
union _dbl_union d;
d.d = val;
if ( (L = word1(d)) !=0)
return P - lo0bits(&L);
L = word0(d) | Exp_msk1;
@ -281,14 +278,14 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
int sudden_underflow;
const char *s, *s0, *s1;
double adj0, tol;
union _dbl_union adj, rv;
Long L;
union _dbl_union adj, rv;
ULong y, z;
Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
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) {
@ -476,20 +473,20 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
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) {
if (e <= Ten_pmax) {
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;
@ -501,17 +498,17 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
*/
e2 = e - i;
e1 -= i;
dval(rv) *= tens[i];
/* rv = */ rounded_product(dval(rv), tens[e2]);
if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
dval(&rv) *= tens[i];
/* rv = */ rounded_product(dval(&rv), tens[e2]);
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;
}
@ -525,48 +522,48 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
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)
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= bigtens[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)
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= tinytens[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];
}
}
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);
@ -763,7 +760,7 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
}
break;
}
if ((dval(adj) = ratio(delta, bs)) <= 2.) {
if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
adj1:
inex = STRTOG_Inexlo;
if (dsign) {
@ -777,15 +774,15 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
irv = STRTOG_Underflow | STRTOG_Inexlo;
break;
}
adj0 = dval(adj) = 1.;
adj0 = dval(&adj) = 1.;
}
else {
adj0 = dval(adj) *= 0.5;
adj0 = dval(&adj) *= 0.5;
if (dsign) {
asub = 0;
inex = STRTOG_Inexlo;
}
if (dval(adj) < 2147483647.) {
if (dval(&adj) < 2147483647.) {
L = adj0;
adj0 -= L;
switch(rd) {
@ -804,12 +801,12 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
inex = STRTOG_Inexact - inex;
}
}
dval(adj) = L;
dval(&adj) = L;
}
}
y = rve + rvbits;
/* adj *= ulp(dval(rv)); */
/* adj *= ulp(&rv); */
/* if (asub) rv -= adj; else rv += adj; */
if (!denorm && rvbits < nbits) {
@ -817,7 +814,7 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
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)
@ -871,15 +868,15 @@ int __strtodg (const char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
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;
}

View File

@ -48,7 +48,7 @@ float __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:

View File

@ -31,33 +31,31 @@ THIS SOFTWARE.
#include "gdtoaimp.h"
double ulp (double x)
double ulp (dbl_union *x)
{
Long L;
union _dbl_union a;
a.d = x;
L = (word0(a) & Exp_mask) - (P-1)*Exp_msk1;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
if (L > 0) {
#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);
}