Aren/NetBSD
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge the new gdtoa: Note that the built-in tests fail the same way as with

Browse Source 
the old gdtoa.
This commit is contained in:
christos 2011-03-20 23:15:35 +00:00
parent f4a9b0b1e6
commit 61e56760bc
37 changed files with 1313 additions and 887 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

201
lib/libc/gdtoa/dtoa.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: dtoa.c,v 1.5 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: dtoa.c,v 1.6 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -68,7 +68,6 @@ THIS SOFTWARE.
*/
#ifdef Honor_FLT_ROUNDS
#define Rounding rounding
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
@ -78,10 +77,10 @@ THIS SOFTWARE.
char *
dtoa
#ifdef KR_headers
(d, mode, ndigits, decpt, sign, rve)
double d; int mode, ndigits, *decpt, *sign; char **rve;
(d0, mode, ndigits, decpt, sign, rve)
double d0; int mode, ndigits, *decpt, *sign; char **rve;
#else
(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
(double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
#endif
{
/* Arguments ndigits, decpt, sign are similar to those
@ -129,14 +128,25 @@ dtoa
#endif
Bigint *b, *b1, *delta, *mhi, *S;
Bigint *mlo = NULL; /* pacify gcc */
double d2, ds, eps;
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) {
@ -144,35 +154,35 @@ dtoa
dtoa_result = 0;
}
#endif
if (word0(d) & Sign_bit) {
d.d = d0;
if (word0(&d) & Sign_bit) {
/* set sign for everything, including 0's and NaNs */
*sign = 1;
word0(d) &= ~Sign_bit; /* clear sign bit */
word0(&d) &= ~Sign_bit; /* clear sign bit */
}
else
*sign = 0;
#if defined(IEEE_Arith) + defined(VAX)
#ifdef IEEE_Arith
if ((word0(d) & Exp_mask) == Exp_mask)
if ((word0(&d) & Exp_mask) == Exp_mask)
#else
if (word0(d) == 0x8000)
if (word0(&d) == 0x8000)
#endif
{
/* Infinity or NaN */
*decpt = 9999;
#ifdef IEEE_Arith
if (!word1(d) && !(word0(d) & 0xfffff))
if (!word1(&d) && !(word0(&d) & 0xfffff))
return nrv_alloc("Infinity", rve, 8);
#endif
return nrv_alloc("NaN", rve, 3);
}
#endif
#ifdef IBM
dval(d) += 0; /* normalize */
dval(&d) += 0; /* normalize */
#endif
if (!dval(d)) {
if (!dval(&d)) {
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
@ -182,37 +192,37 @@ 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);
if (b == NULL)
return NULL;
#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 );
@ -241,21 +251,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;
}
@ -294,10 +304,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;
@ -324,7 +335,7 @@ dtoa
return NULL;
#ifdef Honor_FLT_ROUNDS
if (mode > 1 && rounding != 1)
if (mode > 1 && Rounding != 1)
leftright = 0;
#endif
@ -333,7 +344,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 */
@ -343,7 +354,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 = (unsigned int)j >> 1, i++)
@ -351,32 +362,32 @@ dtoa
ieps++;
ds *= bigtens[i];
}
dval(d) /= ds;
dval(&d) /= ds;
}
else if (( jj1 = -k )!=0) {
dval(d) *= tens[jj1 & 0xf];
dval(&d) *= tens[jj1 & 0xf];
for(j = jj1 >> 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;
}
@ -385,34 +396,34 @@ dtoa
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = 0.5/tens[ilim-1] - dval(eps);
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = dval(d);
dval(d) -= L;
L = dval(&d);
dval(&d) -= L;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps))
if (dval(&d) < dval(&eps))
goto ret1;
if (1. - dval(d) < dval(eps))
if (1. - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d));
if (!(dval(d) -= L))
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d));
if (!(dval(&d) -= L))
ilim = i;
*s++ = '0' + (int)L;
if (i == ilim) {
if (dval(d) > 0.5 + dval(eps))
if (dval(&d) > 0.5 + dval(&eps))
goto bump_up;
else if (dval(d) < 0.5 - dval(eps)) {
else if (dval(&d) < 0.5 - dval(&eps)) {
while(*--s == '0');
s++;
goto ret1;
@ -425,7 +436,7 @@ dtoa
#endif
fast_failed:
s = s0;
dval(d) = dval(d2);
dval(&d) = dval(&d2);
k = k0;
ilim = ilim0;
}
@ -437,22 +448,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
@ -461,13 +472,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) {
@ -544,12 +560,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 */
@ -655,9 +671,9 @@ dtoa
jj1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (jj1 == 0 && mode != 1 && !(word1(d) & 1)
if (jj1 == 0 && mode != 1 && !(word1(&d) & 1)
#ifdef Honor_FLT_ROUNDS
&& rounding >= 1
&& Rounding >= 1
#endif
) {
if (dig == '9')
@ -674,7 +690,7 @@ dtoa
#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) {
@ -685,7 +701,7 @@ dtoa
}
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
switch(rounding) {
switch(Rounding) {
case 0: goto accept_dig;
case 2: goto keep_dig;
}
@ -695,7 +711,11 @@ dtoa
if (b == NULL)
return NULL;
jj1 = cmp(b, S);
#ifdef ROUND_BIASED
if (jjj1 >= 0 /*)*/
#else
if ((jj1 > 0 || (jj1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
}
@ -705,7 +725,7 @@ dtoa
}
if (jj1 > 0) {
#ifdef Honor_FLT_ROUNDS
if (!rounding)
if (!Rounding)
goto accept_dig;
#endif
if (dig == '9') { /* possible if i == 1 */
@ -759,14 +779,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) {
@ -794,9 +819,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)

17
lib/libc/gdtoa/g_Qfmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_Qfmt.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_Qfmt.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -53,15 +53,20 @@ 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 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 };
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;
@ -111,8 +116,8 @@ g_Qfmt(char *buf, void *V, int ndig, unsigned bufsize)
return 0;
mode = 0;
}
s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
return g__fmt(buf, s, se, decpt, sign);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

44
lib/libc/gdtoa/g_ddfmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_ddfmt.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_ddfmt.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,9 +35,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;
@ -45,12 +45,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 < 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]) {
@ -75,7 +91,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)
@ -86,18 +102,18 @@ 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);
z = d2b(dval(&dd[0]), &ex, &bx);
if (z == NULL)
return NULL;
if (dd[1] == 0.)
if (dval(&dd[1]) == 0.)
goto no_y;
x = z;
y = d2b(dd[1], &ey, &by);
y = d2b(dval(&dd[1]), &ey, &by);
if (y == NULL)
return NULL;
if ( (i = ex - ey) !=0) {
@ -159,13 +175,13 @@ 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;
i = STRTOG_Normal;
s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
b = g__fmt(buf, s, se, decpt, z->sign);
b = g__fmt(buf, s, se, decpt, z->sign, bufsize);
if (b == NULL)
return NULL;
Bfree(z);

26
lib/libc/gdtoa/g_dfmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_dfmt.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_dfmt.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,15 +35,20 @@ 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 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 };
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;
@ -54,6 +59,8 @@ g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
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");
}
@ -80,14 +87,13 @@ 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);
if (sign)
i = STRTOG_Normal | STRTOG_Neg;
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
return g__fmt(buf, s, se, decpt, sign);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

17
lib/libc/gdtoa/g_ffmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_ffmt.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_ffmt.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,15 +35,20 @@ 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 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 };
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;
@ -85,8 +90,8 @@ 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);
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
return g__fmt(buf, s, se, decpt, sign);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

17
lib/libc/gdtoa/g_xLfmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_xLfmt.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_xLfmt.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -51,15 +51,20 @@ 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 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 };
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;
@ -105,8 +110,8 @@ g_xLfmt(char *buf, void *V, int ndig, unsigned bufsize)
return 0;
mode = 0;
}
s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
return g__fmt(buf, s, se, decpt, sign);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

17
lib/libc/gdtoa/g_xfmt.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: g_xfmt.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: g_xfmt.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -55,16 +55,21 @@ 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 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 };
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;
@ -111,8 +116,8 @@ g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
return 0;
mode = 0;
}
s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
if (s == NULL)
return NULL;
return g__fmt(buf, s, se, decpt, sign);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

178
lib/libc/gdtoa/gdtoa.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: gdtoa.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: gdtoa.c,v 1.5 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -127,6 +127,7 @@ gdtoa
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 9999.
be = exponent: value = (integer represented by bits) * (2 to the power of be).
mode:
0 ==> shortest string that yields d when read in
@ -161,8 +162,9 @@ gdtoa
int rdir, s2, s5, spec_case, try_quick;
Long L;
Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
double d, d2, ds, eps;
double d2, ds;
char *s, *s0;
U d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
@ -205,21 +207,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 );
@ -239,7 +241,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)
@ -254,13 +256,13 @@ gdtoa
#ifdef IBM
j = be + bbits - 1;
if ( (jj1 = j & 3) !=0)
dval(d) *= 1 << jj1;
word0(d) += j << Exp_shift - 2 & Exp_mask;
dval(&d) *= 1 << jj1;
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;
}
@ -290,11 +292,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;
@ -338,10 +343,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;
@ -352,7 +357,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 /= 2, i++)
@ -364,30 +369,30 @@ gdtoa
else {
ds = 1.;
if ( (jj1 = -k) !=0) {
dval(d) *= tens[jj1 & 0xf];
dval(&d) *= tens[jj1 & 0xf];
for(j = jj1 >> 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;
}
@ -396,42 +401,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;
}
@ -441,7 +444,7 @@ gdtoa
#endif
fast_failed:
s = s0;
dval(d) = d2;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
@ -453,22 +456,22 @@ gdtoa
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = dval(d) / ds;
dval(d) -= L*ds;
for(i = 1;; i++, dval(&d) *= 10.) {
L = dval(&d) / ds;
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (dval(d) == 0.)
if (dval(&d) == 0.)
break;
if (i == ilim) {
if (rdir) {
@ -477,8 +480,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')
@ -489,8 +497,12 @@ gdtoa
}
++*s++;
}
else
else {
inex = STRTOG_Inexlo;
clear_trailing0:
while(*--s == '0'){}
++s;
}
break;
}
}
@ -501,13 +513,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;
@ -583,28 +597,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) {
@ -716,7 +713,11 @@ gdtoa
if (b == NULL)
return NULL;
jj1 = cmp(b, S);
#ifdef ROUND_BIASED
if (jj1 >= 0 /*)*/
#else
if ((jj1 > 0 || (jj1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
@ -780,7 +781,12 @@ gdtoa
if (b == NULL)
return NULL;
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')
@ -796,7 +802,7 @@ gdtoa
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
while(*--s == '0'){}
s++;
++s;
}
ret:
Bfree(S);

22
lib/libc/gdtoa/gdtoa.h
View File

@ -1,4 +1,4 @@
/* $NetBSD: gdtoa.h,v 1.8 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: gdtoa.h,v 1.9 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,6 +35,8 @@ THIS SOFTWARE.
#define GDTOA_H_INCLUDED
#include "arith.h"
#include <stddef.h> /* for size_t */
#include <stdint.h>
#ifndef Long
#define Long int32_t
@ -76,9 +78,9 @@ THIS SOFTWARE.
/* 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,
@ -133,12 +135,12 @@ 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 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_dfmt ANSI((char*, double*, int, size_t));
extern char* g_ffmt ANSI((char*, float*, int, size_t));
extern char* g_Qfmt ANSI((char*, void*, int, size_t));
extern char* g_xfmt ANSI((char*, void*, int, size_t));
extern char* g_xLfmt ANSI((char*, void*, int, size_t));
extern int strtoId ANSI((CONST char*, char**, double*, double*));
extern int strtoIdd ANSI((CONST char*, char**, double*, double*));

88
lib/libc/gdtoa/gdtoaimp.h
View File

@ -1,4 +1,4 @@
/* $NetBSD: gdtoaimp.h,v 1.8 2011/01/21 23:36:49 christos Exp $ */
/* $NetBSD: gdtoaimp.h,v 1.9 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -96,7 +96,12 @@ THIS SOFTWARE.
* #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"
@ -115,7 +120,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,
@ -128,17 +138,22 @@ 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 and spaces;
* if there is only one string of 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.
* NaN(x), where x is a string of hexadecimal digits (optionally
* preceded by 0x or 0X) and spaces; if there is only one string
* of 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.
* Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
* is consumed even when ... has the wrong form (in which case the
* "(...)" is consumed but ignored).
* #define MULTIPLE_THREADS if the system offers preemptively scheduled
* multiple threads. In this case, you must provide (or suitably
* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
@ -150,7 +165,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.
@ -159,11 +174,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.
*/
@ -187,6 +197,9 @@ THIS SOFTWARE.
#define GDTOAIMP_H_INCLUDED
#include "gdtoa.h"
#include "gd_qnan.h"
#ifdef Honor_FLT_ROUNDS
#include <fenv.h>
#endif
#ifdef DEBUG
#include "stdio.h"
@ -281,21 +294,21 @@ typedef union { double d; ULong L[2]; } __attribute__((__may_alias__)) U;
#ifdef YES_ALIAS
#define dval(x) x
#ifdef IEEE_LITTLE_ENDIAN
#define word0(x) ((ULong *)&x)[1]
#define word1(x) ((ULong *)&x)[0]
#define word0(x) ((ULong *)x)[1]
#define word1(x) ((ULong *)x)[0]
#else
#define word0(x) ((ULong *)&x)[0]
#define word1(x) ((ULong *)&x)[1]
#define word0(x) ((ULong *)x)[0]
#define word1(x) ((ULong *)x)[1]
#endif
#else /* !YES_ALIAS */
#ifdef IEEE_LITTLE_ENDIAN
#define word0(x) ( /* LINTED */ (U*)&x)->L[1]
#define word1(x) ( /* LINTED */ (U*)&x)->L[0]
#define word0(x) ( /* LINTED */ (U*)x)->L[1]
#define word1(x) ( /* LINTED */ (U*)x)->L[0]
#else
#define word0(x) ( /* LINTED */ (U*)&x)->L[0]
#define word1(x) ( /* LINTED */ (U*)&x)->L[1]
#define word0(x) ( /* LINTED */ (U*)x)->L[0]
#define word1(x) ( /* LINTED */ (U*)x)->L[1]
#endif
#define dval(x) ( /* LINTED */ (U*)&x)->d
#define dval(x) ( /* LINTED */ (U*)x)->d
#endif /* YES_ALIAS */
/* The following definition of Storeinc is appropriate for MIPS processors.
@ -414,6 +427,11 @@ typedef union { double d; ULong L[2]; } __attribute__((__may_alias__)) 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
@ -563,7 +581,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
extern CONST double bigtens[], tens[], tinytens[];
extern unsigned char hexdig[];
extern Bigint *Balloc ANSI((int));
extern Bigint *Balloc ANSI((size_t));
extern void Bfree ANSI((Bigint*));
extern void ULtof ANSI((ULong*, ULong*, Long, int));
extern void ULtod ANSI((ULong*, ULong*, Long, int));
@ -576,11 +594,11 @@ 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 char *g__fmt ANSI((char*, char*, char*, int, ULong, size_t));
extern int gethex ANSI((CONST char**, CONST FPI*, Long*, Bigint**, int));
extern void hexdig_init_D2A(Void);
extern int hexnan ANSI((CONST char**, CONST FPI*, ULong*));
@ -598,14 +616,16 @@ 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((size_t));
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**, FPI*, Long*, Bigint**, int*));
extern double strtod ANSI((const char *s00, char **se));
extern Bigint *sum ANSI((Bigint*, Bigint*));
extern int trailz ANSI((CONST Bigint*));
extern double ulp ANSI((double));
/*###626 [lint] syntax error '*' [249]%%%*/
/*###626 [lint] incomplete or misplaced function definition [22]%%%*/
extern double ulp ANSI((U*));
#ifdef __cplusplus
}
@ -620,6 +640,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_BIG_ENDIAN
#define _0 0
#define _1 1

167
lib/libc/gdtoa/gethex.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: gethex.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: gethex.c,v 1.5 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -46,20 +46,33 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
#endif
{
Bigint *b;
CONST unsigned char *decpt, *s0, *s, *s1;
int esign, havedig, irv, k, n, nbits, up, zret;
CONST char *decpt, *s, *s0, *s1;
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 char *decimalpoint = localeconv()->decimal_point;
#else
#define decimalpoint '.'
const unsigned char *decimalpoint;
static char *decimalpoint_cache;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = MALLOC(strlen(s0) + 1)) != NULL) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
}
decimalpoint = __UNCONST(s0);
#endif
#endif
if (!hexdig['0'])
if (!hexdig[(unsigned char)'0'])
hexdig_init_D2A();
*bp = 0;
havedig = 0;
s0 = *(CONST unsigned char **)sp + 2;
s0 = *(CONST char **)sp + 2;
while(s0[havedig] == '0')
havedig++;
s0 += havedig;
@ -67,35 +80,54 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
decpt = 0;
zret = 0;
e = 0;
if (!hexdig[*s]) {
if (hexdig[(unsigned char)*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;
if (!hexdig[*s])
#endif
if (!hexdig[(unsigned char)*s])
goto pcheck;
while(*s == '0')
s++;
if (hexdig[*s])
if (hexdig[(unsigned char)*s])
zret = 0;
havedig = 1;
s0 = s;
}
while(hexdig[*s])
while(hexdig[(unsigned char)*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;
while(hexdig[*s])
#endif
while(hexdig[(unsigned char)*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;
@ -103,22 +135,73 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
case '+':
s++;
}
if ((n = hexdig[*s]) == 0 || n > 0x19) {
if ((n = hexdig[(unsigned char)*s]) == 0 || n > 0x19) {
s = s1;
break;
}
e1 = n - 0x10;
while((n = hexdig[*++s]) !=0 && n <= 0x19)
while((n = hexdig[(unsigned char)*++s]) !=0 && n <= 0x19) {
if (e1 & 0xf8000000)
big = 1;
e1 = 10*e1 + n - 0x10;
}
if (esign)
e1 = -e1;
e += e1;
}
*sp = __UNCONST(s);
if (!havedig)
*sp = (char*)__UNCONST(s0) - 1;
if (zret)
return havedig ? STRTOG_Zero : STRTOG_NoNumber;
n = s1 - s0 - 1;
for(k = 0; n > 7; n = (unsigned int)n >> 1)
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 = (unsigned int)nbits >> kshift;
if (nbits & kmask)
++n;
for(j = n, k = 0; (j = (unsigned int)j >> 1) != 0; ++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));
*expt = fpi->emin;
return STRTOG_Normal | STRTOG_Inexlo;
}
n = (int)(s1 - s0) - 1;
for(k = 0; n > (1 << (kshift-2)) - 1; n = (unsigned int)n >> 1)
k++;
b = Balloc(k);
if (b == NULL)
@ -126,20 +209,30 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
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;
}
L |= (hexdig[*s1] & 0x0f) << n;
L |= (hexdig[(unsigned char)*s1] & 0x0f) << n;
n += 4;
}
*x++ = L;
b->wds = n = x - b->x;
n = 32*n - hi0bits(L);
b->wds = n = (int)(x - b->x);
n = ULbits*n - hi0bits(L);
nbits = fpi->nbits;
lostbits = 0;
x = b->x;
@ -150,7 +243,7 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
k = n - 1;
if (x[(unsigned int)k>>kshift] & 1 << (k & kmask)) {
lostbits = 2;
if (k > 1 && any_on(b,k-1))
if (k > 0 && any_on(b,k))
lostbits = 3;
}
}
@ -168,7 +261,10 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, 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;
@ -188,15 +284,22 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
case FPI_Round_down:
if (sign) {
one_bit:
*expt = fpi->emin;
x[0] = b->wds = 1;
dret:
*bp = b;
*expt = 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;
@ -217,7 +320,7 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, 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:
@ -237,7 +340,7 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *expt, Bigint **bp, int sign)
}
else if (b->wds > k
|| ((n = nbits & kmask) !=0
&& hi0bits(x[k-1]) < 32-n)) {
&& hi0bits(x[k-1]) < 32-n)) {
rshift(b,1);
if (++e > fpi->emax)
goto ovfl;

23
lib/libc/gdtoa/hexnan.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: hexnan.c,v 1.3 2006/03/11 18:38:14 kleink Exp $ */
/* $NetBSD: hexnan.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -73,7 +73,13 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
x1 = xe = x;
havedig = hd0 = i = 0;
s = *sp;
while((c = *(CONST unsigned char*)++s) != 0) {
/* allow optional initial 0x or 0X */
while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(CONST unsigned char*)(s+3) > ' ')
s += 2;
while((c = *(CONST unsigned char*)++s)) {
if (!(h = hexdig[c])) {
if (c <= ' ') {
if (hd0 < havedig) {
@ -88,12 +94,25 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
x1 = x;
i = 0;
}
while(*(CONST unsigned char*)(s+1) <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(CONST unsigned char*)(s+3) > ' ')
s += 2;
continue;
}
if (/*(*/ c == ')' && havedig) {
*sp = s + 1;
break;
}
#ifndef GDTOA_NON_PEDANTIC_NANCHECK
do {
if (/*(*/ c == ')') {
*sp = s + 1;
break;
}
} while((c = *++s));
#endif
return STRTOG_NaN;
}
havedig++;

45
lib/libc/gdtoa/makefile
View File

@ -1,4 +1,4 @@
# $NetBSD: makefile,v 1.3 2006/01/25 16:40:57 kleink Exp $
# $NetBSD: makefile,v 1.4 2011/03/20 23:15:35 christos Exp $
# /****************************************************************
# Copyright (C) 1998 by Lucent Technologies
@ -27,11 +27,13 @@
.SUFFIXES: .c .o
CC = cc
CFLAGS = -g -DINFNAN_CHECK
CFLAGS = -g
.c.o:
$(CC) -c $(CFLAGS) $*.c
# invoke "make Printf" to add printf.o to gdtoa.a (if desired)
all: arith.h gd_qnan.h gdtoa.a
arith.h: arithchk.c
@ -44,27 +46,36 @@ gd_qnan.h: arith.h qnan.c
./a.out >gd_qnan.h
rm -f a.out qnan.o
gdtoa.a: dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c g_ffmt.c\
g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c hexnan.c\
misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c\
strtoIx.c strtoIxL.c strtod.c strtodI.c strtodg.c strtof.c strtopQ.c\
strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c\
strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c sum.c ulp.c
gdtoa.a: dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c\
hexnan.c misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c\
strtoIf.c strtoIg.c strtoIx.c strtoIxL.c strtod.c strtodI.c\
strtodg.c strtof.c strtopQ.c strtopd.c strtopdd.c strtopf.c\
strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c strtorf.c\
strtorx.c strtorxL.c sum.c ulp.c
$(CC) -c $(CFLAGS) $?
x=`echo $? | sed 's/\.c/.o/g'` && ar ruv gdtoa.a $$x && rm $$x
ranlib gdtoa.a || true
Printf: all printf.c
$(CC) -c $(CFLAGS) printf.c
ar ruv gdtoa.a printf.o
rm printf.o
touch Printf
# If your system lacks ranlib, you do not need it.
xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h gdtoaimp.h gethex.c\
gmisc.c hd_init.c hexnan.c makefile misc.c qnan.c smisc.c strtoIQ.c\
strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c strtoIxL.c\
strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c strtopQ.c strtopd.c\
strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c\
strtorf.c strtorx.c strtorxL.c sum.c ulp.c
xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c\
g_dfmt.c g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h\
gdtoa_fltrnds.h gdtoaimp.h gethex.c gmisc.c hd_init.c hexnan.c\
makefile misc.c printf.c printf.c0 qnan.c smisc.c stdio1.h\
strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c\
strtoIxL.c strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c\
strtopQ.c strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c\
strtorQ.c strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c\
sum.c ulp.c
# "make xsum.out" to check for transmission errors; source for xsum is
# "make -r xsum.out" to check for transmission errors; source for xsum is
# netlib's "xsum.c from f2c", e.g.,
# ftp://netlib.bell-labs.com/netlib/f2c/xsum.c.gz
@ -73,4 +84,4 @@ xsum.out: xsum0.out $(xs0)
cmp xsum0.out xsum1.out && mv xsum1.out xsum.out || diff xsum[01].out
clean:
rm -f arith.h gd_qnan.h *.[ao] xsum.out xsum1.out
rm -f arith.h gd_qnan.h *.[ao] Printf xsum.out xsum1.out

68
lib/libc/gdtoa/misc.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: misc.c,v 1.5 2009/01/30 23:35:35 lukem Exp $ */
/* $NetBSD: misc.c,v 1.6 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -45,29 +45,31 @@ static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
Bigint *
Balloc
#ifdef KR_headers
(k) int k;
(k) size_t k;
#else
(int k)
(size_t k)
#endif
{
int x;
Bigint *rv;
#ifndef Omit_Private_Memory
unsigned int len;
size_t len;
#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 {
x = 1 << k;
x = 1 << (int)k;
#ifdef Omit_Private_Memory
rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
#else
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
/sizeof(double);
if ((double *)(pmem_next - private_mem + len) <= (double *)PRIVATE_mem) {
if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
rv = (Bigint*)(void *)pmem_next;
pmem_next += len;
}
@ -76,7 +78,7 @@ Balloc
#endif
if (rv == NULL)
return NULL;
rv->k = k;
rv->k = (int)k;
rv->maxwds = x;
}
FREE_DTOA_LOCK(0);
@ -93,10 +95,18 @@ Bfree
#endif
{
if (v) {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
if ((size_t)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);
}
}
}
@ -648,12 +658,12 @@ b2d
{
ULong *xa, *xa0, w, y, z;
int k;
double d;
U d;
#ifdef VAX
ULong d0, d1;
#else
#define d0 word0(d)
#define d1 word1(d)
#define d0 word0(&d)
#define d1 word1(&d)
#endif
xa0 = a->x;
@ -699,10 +709,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
@ -710,12 +720,13 @@ b2d
Bigint *
d2b
#ifdef KR_headers
(d, e, bits) double d; int *e, *bits;
(dd, e, bits) double dd; int *e, *bits;
#else
(double d, int *e, int *bits)
(double dd, int *e, int *bits)
#endif
{
Bigint *b;
U d;
#ifndef Sudden_Underflow
int i;
#endif
@ -723,11 +734,14 @@ d2b
ULong *x, y, z;
#ifdef VAX
ULong d0, d1;
d0 = word0(d) >> 16 | word0(d) << 16;
d1 = word1(d) >> 16 | word1(d) << 16;
#else
#define d0 word0(d)
#define d1 word1(d)
#define d0 word0(&d)
#define d1 word1(&d)
#endif
d.d = dd;
#ifdef VAX
d0 = word0(&d) >> 16 | word0(&d) << 16;
d1 = word1(&d) >> 16 | word1(&d) << 16;
#endif
#ifdef Pack_32
@ -764,10 +778,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
@ -826,7 +836,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;

33
lib/libc/gdtoa/smisc.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: smisc.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: smisc.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -36,9 +36,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;
@ -64,15 +64,15 @@ s2b
i = 9;
if (9 < nd0) {
s += 9;
do {
do {
b = multadd(b, 10, *s++ - '0');
if (b == NULL)
return NULL;
} while(++i < nd0);
s++;
s += dplen;
}
else
s += 10;
s += dplen + 9;
for(; i < nd; i++) {
b = multadd(b, 10, *s++ - '0');
if (b == NULL)
@ -80,7 +80,6 @@ s2b
}
return b;
}
double
ratio
#ifdef KR_headers
@ -89,33 +88,33 @@ ratio
(Bigint *a, Bigint *b)
#endif
{
double da, db;
U da, db;
int k, ka, kb;
dval(da) = b2d(a, &ka);
dval(db) = b2d(b, &kb);
dval(&da) = b2d(a, &ka);
dval(&db) = b2d(b, &kb);
k = ka - kb + ULbits*(a->wds - b->wds);
#ifdef IBM
if (k > 0) {
word0(da) += (k >> 2)*Exp_msk1;
word0(&da) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(da) *= 1 << k;
dval(&da) *= 1 << k;
}
else {
k = -k;
word0(db) += (k >> 2)*Exp_msk1;
word0(&db) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(db) *= 1 << k;
dval(&db) *= 1 << k;
}
#else
if (k > 0)
word0(da) += k*Exp_msk1;
word0(&da) += k*Exp_msk1;
else {
k = -k;
word0(db) += k*Exp_msk1;
word0(&db) += k*Exp_msk1;
}
#endif
return dval(da) / dval(db);
return dval(&da) / dval(&db);
}
#ifdef INFNAN_CHECK

20
lib/libc/gdtoa/strtoIg.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtoIg.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtoIg.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -67,16 +67,20 @@ strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
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);

530
lib/libc/gdtoa/strtod.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtod.c,v 1.5 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtod.c,v 1.6 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -44,16 +44,15 @@ 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
#ifdef Honor_FLT_ROUNDS
#define Rounding rounding
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
@ -65,6 +64,28 @@ __strong_alias(_strtold, strtod)
__weak_alias(strtold, _strtold)
#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
@ -79,20 +100,56 @@ strtod
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
CONST char *s, *s0, *s1;
double aadj, aadj1, adj, rv, rv0;
double aadj;
Long L;
U adj, aadj1, rv, rv0;
ULong y, z;
Bigint *bb = NULL, *bb1, *bd0;
Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
#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;
size_t dplen = strlen(decimalpoint);
#else
char *decimalpoint;
static char *decimalpoint_cache;
static int dplen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
dplen = strlen(s0);
}
decimalpoint = __UNCONST(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 /*}*/
sign = nz0 = nz = decpt = 0;
dval(rv) = 0.;
dval(&rv) = 0.;
for(s = s00;;s++) switch(*s) {
case '-':
sign = 1;
@ -115,7 +172,7 @@ strtod
}
break2:
if (*s == '0') {
#ifndef NO_HEX_FP
#ifndef NO_HEX_FP /*{*/
{
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
Long expt;
@ -124,13 +181,9 @@ strtod
case 'x':
case 'X':
{
#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
#ifdef Honor_FLT_ROUNDS
FPI fpi1 = fpi;
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
@ -151,7 +204,7 @@ strtod
goto ret;
}
}
#endif
#endif /*}*/
nz0 = 1;
while(*++s == '0') ;
if (!*s)
@ -166,13 +219,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++;
@ -201,7 +258,7 @@ strtod
nz = 0;
}
}
}
}/*}*/
dig_done:
e = 0;
if (c == 'e' || c == 'E') {
@ -256,8 +313,8 @@ strtod
--s;
if (!match(&s,"inity"))
++s;
word0(rv) = 0x7ff00000;
word1(rv) = 0;
word0(&rv) = 0x7ff00000;
word1(&rv) = 0;
goto ret;
}
break;
@ -268,13 +325,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
@ -298,13 +355,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
@ -316,6 +373,7 @@ strtod
) {
if (!e)
goto ret;
#ifndef ROUND_BIASED_without_Round_Up
if (e > 0) {
if (e <= Ten_pmax) {
#ifdef VAX
@ -324,11 +382,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
}
@ -340,25 +398,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;
}
@ -368,14 +426,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;
@ -389,12 +448,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*/
@ -403,67 +462,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 = (unsigned int)e1 >> 4;
for(j = 0; e1 > 1; j++, e1 = (unsigned int)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;
@ -472,44 +537,39 @@ strtod
scale = 2*P;
for(j = 0; e1 > 0; j++, e1 = (unsigned int)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.
*/
@ -522,7 +582,7 @@ strtod
/* Put digits into bd: true value = bd * 10^e */
bd0 = s2b(s0, nd0, nd, y);
bd0 = s2b(s0, nd0, nd, y, dplen);
if (bd0 == NULL)
goto ovfl;
@ -531,7 +591,7 @@ strtod
if (bd == NULL)
goto ovfl;
Bcopy(bd, bd0);
bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
if (bb == NULL)
goto ovfl;
bs = i2b(1);
@ -552,16 +612,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
@ -571,7 +638,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
@ -628,7 +695,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) {
@ -638,17 +705,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
@ -657,63 +724,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*/
@ -722,12 +792,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
) {
@ -752,32 +822,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
@ -792,7 +864,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 ==> */
@ -803,10 +875,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
@ -814,16 +886,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
}
@ -835,14 +924,14 @@ strtod
}
if ((aadj = ratio(delta, bs)) <= 2.) {
if (dsign)
aadj = aadj1 = 1.;
else if (word1(rv) || word0(rv) & Bndry_mask) {
aadj = dval(&aadj1) = 1.;
else if (word1(&rv) || word0(&rv) & Bndry_mask) {
#ifndef Sudden_Underflow
if (word1(rv) == Tiny1 && !word0(rv))
if (word1(&rv) == Tiny1 && !word0(&rv))
goto undfl;
#endif
aadj = 1.;
aadj1 = -1.;
dval(&aadj1) = -1.;
}
else {
/* special case -- power of FLT_RADIX to be */
@ -852,45 +941,45 @@ strtod
aadj = 1./FLT_RADIX;
else
aadj *= 0.5;
aadj1 = -aadj;
dval(&aadj1) = -aadj;
}
}
else {
aadj *= 0.5;
aadj1 = dsign ? aadj : -aadj;
dval(&aadj1) = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
switch(Rounding) {
case 2: /* towards +infinity */
aadj1 -= 0.5;
dval(&aadj1) -= 0.5;
break;
case 0: /* towards 0 */
case 3: /* towards -infinity */
aadj1 += 0.5;
dval(&aadj1) += 0.5;
}
#else
if (Flt_Rounds == 0)
aadj1 += 0.5;
dval(&aadj1) += 0.5;
#endif /*Check_FLT_ROUNDS*/
}
y = word0(rv) & Exp_mask;
y = word0(&rv) & Exp_mask;
/* Check for overflow */
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
dval(rv0) = dval(rv);
word0(rv) -= P*Exp_msk1;
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
if ((word0(rv) & Exp_mask) >=
dval(&rv0) = dval(&rv);
word0(&rv) -= P*Exp_msk1;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += dval(&adj);
if ((word0(&rv) & Exp_mask) >=
Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
if (word0(rv0) == Big0 && word1(rv0) == Big1)
if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
goto ovfl;
word0(rv) = Big0;
word1(rv) = Big1;
word0(&rv) = Big0;
word1(&rv) = Big1;
goto cont;
}
else
word0(rv) += P*Exp_msk1;
word0(&rv) += P*Exp_msk1;
}
else {
#ifdef Avoid_Underflow
@ -899,58 +988,58 @@ strtod
if ((z = aadj) == 0)
z = 1;
aadj = z;
aadj1 = dsign ? aadj : -aadj;
dval(&aadj1) = dsign ? aadj : -aadj;
}
word0(aadj1) += (2*P+1)*Exp_msk1 - y;
word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
}
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += dval(&adj);
#else
#ifdef Sudden_Underflow
if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
dval(rv0) = dval(rv);
word0(rv) += P*Exp_msk1;
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
dval(&rv0) = dval(&rv);
word0(&rv) += P*Exp_msk1;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
#ifdef IBM
if ((word0(rv) & Exp_mask) < P*Exp_msk1)
if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
#else
if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
#endif
{
if (word0(rv0) == Tiny0
&& word1(rv0) == Tiny1)
if (word0(&rv0) == Tiny0
&& word1(&rv0) == Tiny1)
goto undfl;
word0(rv) = Tiny0;
word1(rv) = Tiny1;
word0(&rv) = Tiny0;
word1(&rv) = Tiny1;
goto cont;
}
else
word0(rv) -= P*Exp_msk1;
word0(&rv) -= P*Exp_msk1;
}
else {
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
}
#else /*Sudden_Underflow*/
/* Compute adj so that the IEEE rounding rules will
* correctly round rv + adj in some half-way cases.
* If rv * ulp(rv) is denormalized (i.e.,
/* Compute dval(&adj) so that the IEEE rounding rules will
* correctly round rv + dval(&adj) in some half-way cases.
* If rv * ulp(&rv) is denormalized (i.e.,
* y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
* trouble from bits lost to denormalization;
* example: 1.2e-307 .
*/
if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
aadj1 = (double)(int)(aadj + 0.5);
dval(&aadj1) = (double)(int)(aadj + 0.5);
if (!dsign)
aadj1 = -aadj1;
dval(&aadj1) = -dval(&aadj1);
}
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
}
z = word0(rv) & Exp_mask;
z = word0(&rv) & Exp_mask;
#ifndef SET_INEXACT
#ifdef Avoid_Underflow
if (!scale)
@ -960,7 +1049,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;
}
@ -974,12 +1063,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)
@ -987,32 +1081,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 = __UNCONST(s);
return sign ? -dval(rv) : dval(rv);
return sign ? -dval(&rv) : dval(&rv);
}

90
lib/libc/gdtoa/strtodI.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtodI.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtodI.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,16 +35,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;
@ -61,10 +61,6 @@ strtodI(CONST char *s, char **sp, double *dd)
ULong bits[2], sign;
Long exp;
int j, k;
typedef union {
double d[2];
ULong L[4];
} U;
U *u;
k = strtodg(s, sp, &fpi, &exp, bits);
@ -74,17 +70,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
@ -92,80 +88,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;
}

203
lib/libc/gdtoa/strtodg.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtodg.c,v 1.6 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtodg.c,v 1.7 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -93,7 +93,7 @@ increment(Bigint *b)
return b;
}
int
void
#ifdef KR_headers
decrement(b) Bigint *b;
#else
@ -123,7 +123,6 @@ decrement(Bigint *b)
*x++ = y & 0xffff;
} while(borrow && x < xe);
#endif
return STRTOG_Inexlo;
}
static int
@ -179,9 +178,9 @@ set_ones(Bigint *b, int n)
rvOK
#ifdef KR_headers
(d, fpi, expt, bits, exact, rd, irv)
double d; CONST FPI *fpi; Long *expt; ULong *bits; int exact, rd, *irv;
U *d; CONST FPI *fpi; Long *expt; ULong *bits; int exact, rd, *irv;
#else
(double d, CONST FPI *fpi, Long *expt, ULong *bits, int exact, int rd, int *irv)
(U *d, CONST FPI *fpi, Long *expt, ULong *bits, int exact, int rd, int *irv)
#endif
{
Bigint *b;
@ -189,7 +188,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) {
@ -212,9 +211,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;
@ -298,9 +297,9 @@ rvOK
#ifndef VAX
static int
#ifdef KR_headers
mantbits(d) double d;
mantbits(d) U *d;
#else
mantbits(double d)
mantbits(U *d)
#endif
{
ULong L;
@ -335,16 +334,38 @@ strtodg
int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
int sudden_underflow = 0; /* pacify gcc */
CONST char *s, *s0, *s1;
double adj, adj0, rv, tol;
double adj0, tol;
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;
size_t dplen = strlen(decimalpoint);
#else
char *decimalpoint;
static char *decimalpoint_cache;
static int dplen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
dplen = strlen(s0);
}
decimalpoint = __UNCONST(s0);
#endif /*NO_LOCALE_CACHE*/
#else /*USE_LOCALE}{*/
#define dplen 1
#endif /*USE_LOCALE}}*/
e2 = 0; /* XXX gcc */
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) {
@ -399,13 +420,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++;
@ -434,7 +459,7 @@ strtodg
nz = 0;
}
}
}
}/*}*/
dig_done:
e = 0;
if (c == 'e' || c == 'E') {
@ -533,13 +558,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, expt, bits, 1, rd, &irv))
if (rvOK(&rv, fpi, expt, bits, 1, rd, &irv))
goto ret;
}
else if (e > 0) {
@ -547,9 +572,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, expt, bits, i, rd, &irv))
i = fivesbits[e] + mantbits(&rv) <= P;
/* rv = */ rounded_product(dval(&rv), tens[e]);
if (rvOK(&rv, fpi, expt, bits, i, rd, &irv))
goto ret;
e1 -= e;
goto rv_notOK;
@ -562,32 +587,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, expt, bits, 0, rd, &irv))
if (rvOK(&rv, fpi, expt, 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, expt, bits, 0, rd, &irv))
/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
if (rvOK(&rv, fpi, expt, bits, 0, rd, &irv))
goto ret;
e1 -= e;
}
@ -601,45 +626,46 @@ strtodg
e2 = 0;
if (e1 > 0) {
if ( (i = e1 & 15) !=0)
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
if (e1 &= ~15) {
e1 = (unsigned int)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;
for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
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 = (unsigned int)e1 >> 4;
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;
for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
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
@ -650,7 +676,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 */
if (rvb == NULL)
return STRTOG_NoMemory;
rve += e2;
@ -696,7 +722,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);
@ -866,11 +892,8 @@ strtodg
rvb = increment(rvb);
if (rvb == NULL)
return STRTOG_NoMemory;
if ( (j = rvbits & kmask) !=0)
j = ULbits - j;
if (hi0bits(rvb->x[(unsigned int)(rvb->wds - 1)
>> kshift])
!= j)
j = kmask & (ULbits - (rvbits & kmask));
if (hi0bits(rvb->x[rvb->wds - 1]) != j)
rvbits++;
irv = STRTOG_Normal | STRTOG_Inexhi;
}
@ -882,7 +905,7 @@ strtodg
}
break;
}
if ((dval(adj) = ratio(delta, bs)) <= 2.) {
if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
adj1:
inex = STRTOG_Inexlo;
if (dsign) {
@ -896,15 +919,15 @@ strtodg
irv = STRTOG_Underflow | STRTOG_Inexlo;
break;
}
adj0 = dval(adj) = 1.;
adj0 = dval(&adj) = 1.;
}
else {
adj0 = dval(adj) *= 0.5;
adj0 = dval(&adj) *= 0.5;
if (dsign) {
asub = 0;
inex = STRTOG_Inexlo;
}
if (dval(adj) < 2147483647.) {
if (dval(&adj) < 2147483647.) {
L = adj0;
adj0 -= L;
switch(rd) {
@ -923,12 +946,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) {
@ -938,7 +961,7 @@ strtodg
rve -= j;
rvbits = nbits;
}
ab = d2b(dval(adj), &abe, &abits);
ab = d2b(dval(&adj), &abe, &abits);
if (ab == NULL)
return STRTOG_NoMemory;
if (abe < 0)
@ -1000,15 +1023,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;
}
@ -1033,6 +1056,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;
*expt = 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;
@ -1049,12 +1095,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)

10
lib/libc/gdtoa/strtodnrp.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtodnrp.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtodnrp.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -46,13 +46,13 @@ strtod(s, sp) CONST char *s; char **sp;
strtod(CONST char *s, char **sp)
#endif
{
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
static const FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
ULong bits[2];
Long exp;
Long expt;
int k;
union { ULong L[2]; double d; } u;
k = strtodg(s, sp, &fpi, &exp, bits);
k = strtodg(s, sp, &fpi, &expt, bits);
if (k == STRTOG_NoMemory) {
errno = ERANGE;
u.L[0] = Big0;
@ -67,7 +67,7 @@ strtod(CONST char *s, char **sp)
case STRTOG_Normal:
u.L[_1] = bits[0];
u.L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
u.L[_0] = (bits[1] & ~0x100000) | ((expt + 0x3ff + 52) << 20);
break;
case STRTOG_Denormal:

11
lib/libc/gdtoa/strtof.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtof.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtof.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -45,13 +45,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 CONST FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
ULong bits[1];
Long expt;
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, &expt, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory) {
errno = ERANGE;
return HUGE_VALF;

11
lib/libc/gdtoa/strtopQ.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopQ.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopQ.c,v 1.5 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -58,13 +58,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 CONST FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI };
ULong bits[4];
Long expt;
int k;
ULong *L = (ULong*)V;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
k = strtodg(s, sp, &fpi, &expt, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
switch(k & STRTOG_Retmask) {

13
lib/libc/gdtoa/strtopd.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopd.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopd.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -42,12 +42,17 @@ strtopd(CONST char *s, char **sp, double *d)
{
static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
ULong bits[2];
Long exp;
Long expt;
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, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtod((ULong*)d, bits, exp, k);
ULtod((ULong*)d, bits, expt, k);
return k;
}

69
lib/libc/gdtoa/strtopdd.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopdd.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopdd.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -41,20 +41,25 @@ strtopdd(CONST char *s, char **sp, double *dd)
#endif
{
#ifdef Sudden_Underflow
static 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 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;
Long expt;
int i, j, rv;
typedef union {
double d[2];
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, &expt, bits);
if (rv == STRTOG_NoMemory)
return rv;
u = (U*)dd;
@ -66,36 +71,36 @@ 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;
exp += 0x3ff + 52;
u->L[_0] = (bits[2] >> 21) | ((bits[3] << 11) & 0xfffff)
| ((expt + 0x3ff + 105) << 20);
expt += 0x3ff + 52;
if (bits[1] &= 0x1fffff) {
i = hi0bits(bits[1]) - 11;
if (i >= exp) {
i = exp - 1;
exp = 0;
if (i >= expt) {
i = expt - 1;
expt = 0;
}
else
exp -= i;
expt -= 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;
}
}
else if (bits[0]) {
i = hi0bits(bits[0]) + 21;
if (i >= exp) {
i = exp - 1;
exp = 0;
if (i >= expt) {
i = expt - 1;
expt = 0;
}
else
exp -= i;
expt -= 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;
}
}
@ -104,7 +109,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) | (expt << 20);
break;
case STRTOG_Denormal:
@ -123,10 +128,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;
@ -135,34 +140,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:

15
lib/libc/gdtoa/strtopf.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopf.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopf.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -40,12 +40,17 @@ strtopf(s, sp, f) CONST char *s; char **sp; float *f;
strtopf(CONST char *s, char **sp, float *f)
#endif
{
static 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;
Long expt;
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, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
L = (ULong*)f;
@ -57,7 +62,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) | ((expt + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:

14
lib/libc/gdtoa/strtopx.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopx.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopx.c,v 1.5 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -60,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 const FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
ULong bits[2];
Long expt;
int k;
UShort *L = (UShort*)V;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
k = strtodg(s, sp, &fpi, &expt, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
switch(k & STRTOG_Retmask) {
@ -91,7 +96,8 @@ strtopx(CONST char *s, char **sp, void *V)
case STRTOG_Infinite:
L[_0] = 0x7fff;
L[_1] = L[_2] = L[_3] = L[_4] = 0;
L[_1] = 0x8000;
L[_2] = L[_3] = L[_4] = 0;
break;
case STRTOG_NaN:

14
lib/libc/gdtoa/strtopxL.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtopxL.c,v 1.4 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtopxL.c,v 1.5 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -56,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 CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
ULong bits[2];
Long expt;
int k;
ULong *L = (ULong*)V;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
k = strtodg(s, sp, &fpi, &expt, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
switch(k & STRTOG_Retmask) {
@ -81,7 +86,8 @@ strtopxL(CONST char *s, char **sp, void *V)
case STRTOG_Infinite:
L[_0] = 0x7fff << 16;
L[_1] = L[_2] = 0;
L[_1] = 0x80000000;
L[_2] = 0;
break;
case STRTOG_NaN:

14
lib/libc/gdtoa/strtorQ.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtorQ.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtorQ.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -53,9 +53,9 @@ THIS SOFTWARE.
void
#ifdef KR_headers
ULtoQ(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
ULtoQ(L, bits, expt, k) ULong *L; ULong *bits; Long expt; int k;
#else
ULtoQ(ULong *L, ULong *bits, Long exp, int k)
ULtoQ(ULong *L, ULong *bits, Long expt, int k)
#endif
{
switch(k & STRTOG_Retmask) {
@ -69,7 +69,7 @@ ULtoQ(ULong *L, ULong *bits, Long exp, int k)
L[_3] = bits[0];
L[_2] = bits[1];
L[_1] = bits[2];
L[_0] = (bits[3] & ~0x10000) | ((exp + 0x3fff + 112) << 16);
L[_0] = (bits[3] & ~0x10000) | ((expt + 0x3fff + 112) << 16);
break;
case STRTOG_Denormal:
@ -104,7 +104,7 @@ strtorQ(CONST char *s, char **sp, int rounding, void *L)
static FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
FPI *fpi, fpi1;
ULong bits[4];
Long exp;
Long expt;
int k;
fpi = &fpi0;
@ -113,9 +113,9 @@ strtorQ(CONST char *s, char **sp, int rounding, void *L)
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtoQ((ULong*)L, bits, exp, k);
ULtoQ((ULong*)L, bits, expt, k);
return k;
}

71
lib/libc/gdtoa/strtordd.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtordd.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtordd.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,9 +35,9 @@ THIS SOFTWARE.
void
#ifdef KR_headers
ULtodd(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
ULtodd(L, bits, expt, k) ULong *L; ULong *bits; Long expt; int k;
#else
ULtodd(ULong *L, ULong *bits, Long exp, int k)
ULtodd(ULong *L, ULong *bits, Long expt, int k)
#endif
{
int i, j;
@ -50,36 +50,36 @@ 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;
exp += 0x3ff + 52;
L[_0] = (bits[2] >> 21) | (bits[3] << 11 & 0xfffff)
| ((expt + 0x3ff + 105) << 20);
expt += 0x3ff + 52;
if (bits[1] &= 0x1fffff) {
i = hi0bits(bits[1]) - 11;
if (i >= exp) {
i = exp - 1;
exp = 0;
if (i >= expt) {
i = expt - 1;
expt = 0;
}
else
exp -= i;
expt -= 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;
}
}
else if (bits[0]) {
i = hi0bits(bits[0]) + 21;
if (i >= exp) {
i = exp - 1;
exp = 0;
if (i >= expt) {
i = expt - 1;
expt = 0;
}
else
exp -= i;
expt -= 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;
}
}
@ -88,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) | (expt << 20);
break;
case STRTOG_Denormal:
@ -107,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;
@ -119,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:
@ -180,13 +180,14 @@ strtordd(CONST char *s, char **sp, int rounding, double *dd)
#endif
{
#ifdef Sudden_Underflow
static FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
static CONST FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
#else
static FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
static CONST FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
#endif
FPI *fpi, fpi1;
CONST FPI *fpi;
FPI fpi1;
ULong bits[4];
Long exp;
Long expt;
int k;
fpi = &fpi0;
@ -195,9 +196,9 @@ strtordd(CONST char *s, char **sp, int rounding, double *dd)
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtodd((ULong*)dd, bits, exp, k);
ULtodd((ULong*)dd, bits, expt, k);
return k;
}

19
lib/libc/gdtoa/strtorf.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtorf.c,v 1.2 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtorf.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -35,9 +35,9 @@ THIS SOFTWARE.
void
#ifdef KR_headers
ULtof(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
ULtof(L, bits, expt, k) ULong *L; ULong *bits; Long expt; int k;
#else
ULtof(ULong *L, ULong *bits, Long exp, int k)
ULtof(ULong *L, ULong *bits, Long expt, int k)
#endif
{
switch(k & STRTOG_Retmask) {
@ -48,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) | ((expt + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:
@ -73,10 +73,11 @@ 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 FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
FPI *fpi, fpi1;
static CONST FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
CONST FPI *fpi;
FPI fpi1;
ULong bits[1];
Long exp;
Long expt;
int k;
fpi = &fpi0;
@ -85,9 +86,9 @@ strtorf(CONST char *s, char **sp, int rounding, float *f)
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtof((ULong*)f, bits, exp, k);
ULtof((ULong*)f, bits, expt, k);
return k;
}

22
lib/libc/gdtoa/strtorx.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtorx.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtorx.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -55,9 +55,9 @@ THIS SOFTWARE.
void
#ifdef KR_headers
ULtox(L, bits, exp, k) UShort *L; ULong *bits; Long exp; int k;
ULtox(L, bits, expt, k) UShort *L; ULong *bits; Long expt; int k;
#else
ULtox(UShort *L, ULong *bits, Long exp, int k)
ULtox(UShort *L, ULong *bits, Long expt, int k)
#endif
{
switch(k & STRTOG_Retmask) {
@ -72,7 +72,7 @@ ULtox(UShort *L, ULong *bits, Long exp, int k)
case STRTOG_Normal:
case STRTOG_NaNbits:
L[_0] = exp + 0x3fff + 63;
L[_0] = expt + 0x3fff + 63;
normal_bits:
L[_4] = (UShort)bits[0];
L[_3] = (UShort)(bits[0] >> 16);
@ -82,7 +82,8 @@ ULtox(UShort *L, ULong *bits, Long exp, int k)
case STRTOG_Infinite:
L[_0] = 0x7fff;
L[_1] = L[_2] = L[_3] = L[_4] = 0;
L[_1] = 0x8000;
L[_2] = L[_3] = L[_4] = 0;
break;
case STRTOG_NaN:
@ -103,10 +104,11 @@ 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 FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
FPI *fpi, fpi1;
static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
CONST FPI *fpi;
FPI fpi1;
ULong bits[2];
Long exp;
Long expt;
int k;
fpi = &fpi0;
@ -115,9 +117,9 @@ strtorx(CONST char *s, char **sp, int rounding, void *L)
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtox((UShort*)L, bits, exp, k);
ULtox((UShort*)L, bits, expt, k);
return k;
}

22
lib/libc/gdtoa/strtorxL.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: strtorxL.c,v 1.3 2008/03/21 23:13:48 christos Exp $ */
/* $NetBSD: strtorxL.c,v 1.4 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -51,9 +51,9 @@ THIS SOFTWARE.
void
#ifdef KR_headers
ULtoxL(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
ULtoxL(L, bits, expt, k) ULong *L; ULong *bits; Long expt; int k;
#else
ULtoxL(ULong *L, ULong *bits, Long exp, int k)
ULtoxL(ULong *L, ULong *bits, Long expt, int k)
#endif
{
switch(k & STRTOG_Retmask) {
@ -65,14 +65,15 @@ ULtoxL(ULong *L, ULong *bits, Long exp, int k)
case STRTOG_Normal:
case STRTOG_Denormal:
case STRTOG_NaNbits:
L[_0] = (exp + 0x3fff + 63) << 16;
L[_0] = (expt + 0x3fff + 63) << 16;
L[_1] = bits[1];
L[_2] = bits[0];
break;
case STRTOG_Infinite:
L[_0] = 0x7fff << 16;
L[_1] = L[_2] = 0;
L[_1] = 0x80000000;
L[_2] = 0;
break;
case STRTOG_NaN:
@ -91,10 +92,11 @@ 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 FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
FPI *fpi, fpi1;
static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
CONST FPI *fpi;
FPI fpi1;
ULong bits[2];
Long exp;
Long expt;
int k;
fpi = &fpi0;
@ -103,9 +105,9 @@ strtorxL(CONST char *s, char **sp, int rounding, void *L)
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
k = strtodg(s, sp, fpi, &expt, bits);
if (k == STRTOG_NoMemory)
return k;
ULtoxL((ULong*)L, bits, exp, k);
ULtoxL((ULong*)L, bits, expt, k);
return k;
}

6
lib/libc/gdtoa/test/Qtest.c
View File

@ -60,15 +60,15 @@ THIS SOFTWARE.
#undef _0
#undef _1
/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
/* one or the other of IEEE_BIG_ENDIAN or IEEE_LITTLE_ENDIAN should be #defined */
#ifdef IEEE_MC68k
#ifdef IEEE_BIG_ENDIAN
#define _0 0
#define _1 1
#define _2 2
#define _3 3
#endif
#ifdef IEEE_8087
#ifdef IEEE_LITTLE_ENDIAN
#define _0 3
#define _1 2
#define _2 1

49
lib/libc/gdtoa/test/makefile
View File

@ -23,17 +23,15 @@
#
# ****************************************************************/
.SUFFIXES: .c .o
CC = cc
CFLAGS = -g -I..
A = ../gdtoa.a
CFLAGS+=-I${.CURDIR}/../../arch/${MACHINE_ARCH}/gdtoa
L = -lm
L1 =
L1 = -lm
#use "L1=-lm" when compiled with -DHonor_FLTP_ROUNDS or -DUSE_LOCALE
INFFIX = | sed 's/[Ii][Nn][Ff][intyINTY]*/Infinity/g'
.c.o:
$(CC) -c $(CFLAGS) $*.c
.PATH: /usr/src/lib/libc/gdtoa
all: dt dItest ddtest dtest ftest Qtest xLtest xtest ddtestsi dItestsi tests
@ -41,49 +39,48 @@ dt = dt.o $A
dt: $(dt)
$(CC) -o dt $(dt) $L
dItest = dItest.o getround.o $A
dItest = dItest.o getround.o g_dfmt.o strtoId.o strtodI.o \
g__fmt.o strtoIg.o
dItest: $(dItest)
$(CC) -o dItest $(dItest) $(L1)
ddtest = ddtest.o getround.o $A
ddtest = ddtest.o getround.o g_dfmt.o strtordd.o strtopdd.o g_ddfmt.o \
strtoIdd.o g__fmt.o strtoIg.o
ddtest: $(ddtest)
$(CC) -o ddtest $(ddtest) $L
dtest = dtest.o getround.o $A
dtest = dtest.o getround.o g_dfmt.o strtopd.o strtoId.o \
g__fmt.o strtoIg.o
dtest: $(dtest)
$(CC) -o dtest $(dtest) $L
ftest = ftest.o getround.o $A
ftest = ftest.o getround.o strtorf.o strtopf.o g_ffmt.o strtoIf.o \
g__fmt.o strtoIg.o
ftest: $(ftest)
$(CC) -o ftest $(ftest) $(L1)
Qtest = Qtest.o getround.o $A
Qtest = Qtest.o getround.o strtorQ.o g_Qfmt.o strtoIQ.o strtopQ.o \
g__fmt.o strtoIg.o
Qtest: $(Qtest)
$(CC) -o Qtest $(Qtest) $(L1)
xtest = xtest.o getround.o $A
xtest = xtest.o getround.o strtorx.o g_xfmt.o strtoIx.o \
g__fmt.o strtoIg.o
xtest: $(xtest)
$(CC) -o xtest $(xtest) $(L1)
xLtest = xLtest.o getround.o $A
xLtest = xLtest.o getround.o strtorxL.o strtopxL.o g_xLfmt.o strtoIxL.o \
g__fmt.o strtoIg.o
xLtest: $(xLtest)
$(CC) -o xLtest $(xLtest) $(L1)
strtopddSI.o: strtopddSI.c ../strtopdd.c
strtorddSI.o: strtorddSI.c ../strtordd.c
strtodISI.o: strtodISI.c ../strtodI.c
strtoIddSI.o: strtoIddSI.c ../strtoIdd.c
strtoIdSI.o: strtoIdSI.c ../strtoId.c
ddtestsi = ddtest.o strtopddSI.o strtorddSI.o strtoIddSI.o getround.o $A
ddtestsi = ddtest.o strtopddSI.o strtorddSI.o strtoIddSI.o getround.o \
g_dfmt.o g_ddfmt.o g__fmt.o strtoIg.o
ddtestsi: $(ddtestsi)
$(CC) -o ddtestsi $(ddtestsi) $L
dItestsi = dItest.o strtodISI.o strtoIdSI.o getround.o $A
dItestsi = dItest.o strtodISI.o strtoIdSI.o getround.o \
g_dfmt.o g__fmt.o strtoIg.o
dItestsi: $(dItestsi)
$(CC) -o dItestsi $(dItestsi) $(L1)
@ -92,7 +89,7 @@ strtodt: $(strtodt)
$(CC) -o strtodt $(strtodt) $L
pftest = pftest.o $A
pftest: ../Printf $(pftest)
pftest: $(pftest)
$(CC) -o pftest $(pftest) $L
## On Intel (and Intel-like) systems using extended-precision registers
@ -102,7 +99,7 @@ pftest: ../Printf $(pftest)
## (which does all computations in integer arithmetic) and should show no
## unexpected results.
strtodtnrp = strtodt.o ../strtodnrp.c $A
strtodtnrp = strtodt.o strtodnrp.o
strtodtnrp: $(strtodtnrp)
$(CC) -o strtodtnrp $(strtodtnrp)

6
lib/libc/gdtoa/test/xLtest.c
View File

@ -63,14 +63,14 @@ THIS SOFTWARE.
#undef _0
#undef _1
/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
/* one or the other of IEEE_BIG_ENDIAN or IEEE_LITTLE_ENDIAN should be #defined */
#ifdef IEEE_MC68k
#ifdef IEEE_BIG_ENDIAN
#define _0 0
#define _1 1
#define _2 2
#endif
#ifdef IEEE_8087
#ifdef IEEE_LITTLE_ENDIAN
#define _0 2
#define _1 1
#define _2 0

6
lib/libc/gdtoa/test/xtest.c
View File

@ -61,16 +61,16 @@ THIS SOFTWARE.
#undef _0
#undef _1
/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
/* one or the other of IEEE_BIG_ENDIAN or IEEE_LITTLE_ENDIAN should be #defined */
#ifdef IEEE_MC68k
#ifdef IEEE_BIG_ENDIAN
#define _0 0
#define _1 1
#define _2 2
#define _3 3
#define _4 4
#endif
#ifdef IEEE_8087
#ifdef IEEE_LITTLE_ENDIAN
#define _0 4
#define _1 3
#define _2 2

22
lib/libc/gdtoa/ulp.c
View File

@ -1,4 +1,4 @@
/* $NetBSD: ulp.c,v 1.2 2006/01/25 15:27:42 kleink Exp $ */
/* $NetBSD: ulp.c,v 1.3 2011/03/20 23:15:35 christos Exp $ */
/****************************************************************
@ -36,13 +36,13 @@ THIS SOFTWARE.
double
ulp
#ifdef KR_headers
(x) double x;
(x) U *x;
#else
(double x)
(U *x)
#endif
{
Long L;
double a;
U a;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
@ -51,22 +51,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 = (unsigned int)-L >> Exp_shift;
if (L < Exp_shift) {
word0(a) = 0x80000 >> L;
word1(a) = 0;
word0(&a) = 0x80000 >> L;
word1(&a) = 0;
}
else {
word0(a) = 0;
word0(&a) = 0;
L -= Exp_shift;
word1(a) = L >= 31 ? 1 : 1 << (31 - L);
word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
}
}
#endif
return a;
return dval(&a);
}