ref: 403fef45c1a65aa74c449764e2857597fc03290a
dir: /sys/src/ape/lib/ap/stdio/strtod.c/
#include "fconv.h"
/* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg).
*
* This strtod returns a nearest machine number to the input decimal
* string (or sets errno to ERANGE). With IEEE arithmetic, ties are
* broken by the IEEE round-even rule. Otherwise ties are broken by
* biased rounding (add half and chop).
*
* Inspired loosely by William D. Clinger's paper "How to Read Floating
* Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
*
* Modifications:
*
* 1. We only require IEEE, IBM, or VAX double-precision
* arithmetic (not IEEE double-extended).
* 2. We get by with floating-point arithmetic in a case that
* Clinger missed -- when we're computing d * 10^n
* for a small integer d and the integer n is not too
* much larger than 22 (the maximum integer k for which
* we can represent 10^k exactly), we may be able to
* compute (d*10^k) * 10^(e-k) with just one roundoff.
* 3. Rather than a bit-at-a-time adjustment of the binary
* result in the hard case, we use floating-point
* arithmetic to determine the adjustment to within
* one bit; only in really hard cases do we need to
* compute a second residual.
* 4. Because of 3., we don't need a large table of powers of 10
* for ten-to-e (just some small tables, e.g. of 10^k
* for 0 <= k <= 22).
*/
#ifdef RND_PRODQUOT
#define rounded_product(a,b) a = rnd_prod(a, b)
#define rounded_quotient(a,b) a = rnd_quot(a, b)
extern double rnd_prod(double, double), rnd_quot(double, double);
#else
#define rounded_product(a,b) a *= b
#define rounded_quotient(a,b) a /= b
#endif
static double
ulp(double xarg)
{
register long L;
Dul a;
Dul x;
x.d = xarg;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
if (L > 0) {
#endif
#ifdef IBM
L |= Exp_msk1 >> 4;
#endif
word0(a) = L;
word1(a) = 0;
#ifndef Sudden_Underflow
}
else {
L = -L >> Exp_shift;
if (L < Exp_shift) {
word0(a) = 0x80000 >> L;
word1(a) = 0;
}
else {
word0(a) = 0;
L -= Exp_shift;
word1(a) = L >= 31 ? 1 : 1 << 31 - L;
}
}
#endif
return a.d;
}
static Bigint *
s2b(CONST char *s, int nd0, int nd, unsigned long y9)
{
Bigint *b;
int i, k;
long x, y;
x = (nd + 8) / 9;
for(k = 0, y = 1; x > y; y <<= 1, k++) ;
#ifdef Pack_32
b = Balloc(k);
b->x[0] = y9;
b->wds = 1;
#else
b = Balloc(k+1);
b->x[0] = y9 & 0xffff;
b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
#endif
i = 9;
if (9 < nd0) {
s += 9;
do b = multadd(b, 10, *s++ - '0');
while(++i < nd0);
s++;
}
else
s += 10;
for(; i < nd; i++)
b = multadd(b, 10, *s++ - '0');
return b;
}
static double
b2d(Bigint *a, int *e)
{
unsigned long *xa, *xa0, w, y, z;
int k;
Dul d;
#ifdef VAX
unsigned long d0, d1;
#else
#define d0 word0(d)
#define d1 word1(d)
#endif
xa0 = a->x;
xa = xa0 + a->wds;
y = *--xa;
#ifdef DEBUG
if (!y) Bug("zero y in b2d");
#endif
k = hi0bits(y);
*e = 32 - k;
#ifdef Pack_32
if (k < Ebits) {
d0 = Exp_1 | y >> Ebits - k;
w = xa > xa0 ? *--xa : 0;
d1 = y << (32-Ebits) + k | w >> Ebits - k;
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
if (k -= Ebits) {
d0 = Exp_1 | y << k | z >> 32 - k;
y = xa > xa0 ? *--xa : 0;
d1 = z << k | y >> 32 - k;
}
else {
d0 = Exp_1 | y;
d1 = z;
}
#else
if (k < Ebits + 16) {
z = xa > xa0 ? *--xa : 0;
d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
w = xa > xa0 ? *--xa : 0;
y = xa > xa0 ? *--xa : 0;
d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
w = xa > xa0 ? *--xa : 0;
k -= Ebits + 16;
d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
y = xa > xa0 ? *--xa : 0;
d1 = w << k + 16 | y << k;
#endif
ret_d:
#ifdef VAX
word0(d) = d0 >> 16 | d0 << 16;
word1(d) = d1 >> 16 | d1 << 16;
#else
#undef d0
#undef d1
#endif
return d.d;
}
static double
ratio(Bigint *a, Bigint *b)
{
Dul da, db;
int k, ka, kb;
da.d = b2d(a, &ka);
db.d = b2d(b, &kb);
#ifdef Pack_32
k = ka - kb + 32*(a->wds - b->wds);
#else
k = ka - kb + 16*(a->wds - b->wds);
#endif
#ifdef IBM
if (k > 0) {
word0(da) += (k >> 2)*Exp_msk1;
if (k &= 3)
da *= 1 << k;
}
else {
k = -k;
word0(db) += (k >> 2)*Exp_msk1;
if (k &= 3)
db *= 1 << k;
}
#else
if (k > 0)
word0(da) += k*Exp_msk1;
else {
k = -k;
word0(db) += k*Exp_msk1;
}
#endif
return da.d / db.d;
}
double
strtod(CONST char *s00, char **se)
{
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
CONST char *s, *s0, *s1;
double aadj, aadj1, adj;
Dul rv, rv0;
long L;
unsigned long y, z;
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
sign = nz0 = nz = 0;
rv.d = 0.;
for(s = s00;;s++) switch(*s) {
case '-':
sign = 1;
/* no break */
case '+':
if (*++s)
goto break2;
/* no break */
case 0:
s = s00;
goto ret;
case '\t':
case '\n':
case '\v':
case '\f':
case '\r':
case ' ':
continue;
default:
goto break2;
}
break2:
if (*s == '0') {
nz0 = 1;
while(*++s == '0') ;
if (!*s)
goto ret;
}
s0 = s;
y = z = 0;
for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
if (nd < 9)
y = 10*y + c - '0';
else if (nd < 16)
z = 10*z + c - '0';
nd0 = nd;
if (c == '.') {
c = *++s;
if (!nd) {
for(; c == '0'; c = *++s)
nz++;
if (c > '0' && c <= '9') {
s0 = s;
nf += nz;
nz = 0;
goto have_dig;
}
goto dig_done;
}
for(; c >= '0' && c <= '9'; c = *++s) {
have_dig:
nz++;
if (c -= '0') {
nf += nz;
for(i = 1; i < nz; i++)
if (nd++ < 9)
y *= 10;
else if (nd <= DBL_DIG + 1)
z *= 10;
if (nd++ < 9)
y = 10*y + c;
else if (nd <= DBL_DIG + 1)
z = 10*z + c;
nz = 0;
}
}
}
dig_done:
e = 0;
if (c == 'e' || c == 'E') {
if (!nd && !nz && !nz0) {
s = s00;
goto ret;
}
s00 = s;
esign = 0;
switch(c = *++s) {
case '-':
esign = 1;
case '+':
c = *++s;
}
if (c >= '0' && c <= '9') {
while(c == '0')
c = *++s;
if (c > '0' && c <= '9') {
e = c - '0';
s1 = s;
while((c = *++s) >= '0' && c <= '9')
e = 10*e + c - '0';
if (s - s1 > 8)
/* Avoid confusion from exponents
* so large that e might overflow.
*/
e = 9999999;
if (esign)
e = -e;
}
else
e = 0;
}
else
s = s00;
}
if (!nd) {
if (!nz && !nz0)
s = s00;
goto ret;
}
e1 = e -= nf;
/* Now we have nd0 digits, starting at s0, followed by a
* decimal point, followed by nd-nd0 digits. The number we're
* after is the integer represented by those digits times
* 10**e */
if (!nd0)
nd0 = nd;
k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
rv.d = y;
if (k > 9)
rv.d = tens[k - 9] * rv.d + z;
bd0 = 0;
if (nd <= DBL_DIG
#ifndef RND_PRODQUOT
&& FLT_ROUNDS == 1
#endif
) {
if (!e)
goto ret;
if (e > 0) {
if (e <= Ten_pmax) {
#ifdef VAX
goto vax_ovfl_check;
#else
/* rv = */ rounded_product(rv.d, tens[e]);
goto ret;
#endif
}
i = DBL_DIG - nd;
if (e <= Ten_pmax + i) {
/* A fancier test would sometimes let us do
* this for larger i values.
*/
e -= i;
rv.d *= 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(rv.d, tens[e]);
if ((word0(rv) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
goto ovfl;
word0(rv) += P*Exp_msk1;
#else
/* rv = */ rounded_product(rv.d, tens[e]);
#endif
goto ret;
}
}
else if (e >= -Ten_pmax) {
/* rv = */ rounded_quotient(rv.d, tens[-e]);
goto ret;
}
}
e1 += nd - k;
/* Get starting approximation = rv * 10**e1 */
if (e1 > 0) {
if (nd0 + e1 - 1 > DBL_MAX_10_EXP)
goto ovfl;
if (i = e1 & 15)
rv.d *= tens[i];
if (e1 &= ~15) {
if (e1 > DBL_MAX_10_EXP) {
ovfl:
errno = ERANGE;
rv.d = HUGE_VAL;
if (bd0)
goto retfree;
goto ret;
}
if (e1 >>= 4) {
for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1)
rv.d *= bigtens[j];
/* The last multiplication could overflow. */
word0(rv) -= P*Exp_msk1;
rv.d *= 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;
}
else
word0(rv) += P*Exp_msk1;
}
}
}
else if (e1 < 0) {
e1 = -e1;
if (i = e1 & 15)
rv.d /= tens[i];
if (e1 &= ~15) {
e1 >>= 4;
if (e1 >= 1 << n_bigtens)
goto undfl;
for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1)
rv.d *= tinytens[j];
/* The last multiplication could underflow. */
rv0.d = rv.d;
rv.d *= tinytens[j];
if (rv.d == 0) {
rv.d = 2.*rv0.d;
rv.d *= tinytens[j];
if (rv.d == 0) {
undfl:
rv.d = 0.;
errno = ERANGE;
if (bd0)
goto retfree;
goto ret;
}
word0(rv) = Tiny0;
word1(rv) = Tiny1;
/* The refinement below will clean
* this approximation up.
*/
}
}
}
/* Now the hard part -- adjusting rv to the correct value.*/
/* Put digits into bd: true value = bd * 10^e */
bd0 = s2b(s0, nd0, nd, y);
for(;;) {
bd = Balloc(bd0->k);
Bcopy(bd, bd0);
bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */
bs = i2b(1);
if (e >= 0) {
bb2 = bb5 = 0;
bd2 = bd5 = e;
}
else {
bb2 = bb5 = -e;
bd2 = bd5 = 0;
}
if (bbe >= 0)
bb2 += bbe;
else
bd2 -= bbe;
bs2 = bb2;
#ifdef Sudden_Underflow
#ifdef IBM
j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
#else
j = P + 1 - bbbits;
#endif
#else
i = bbe + bbbits - 1; /* logb(rv) */
if (i < Emin) /* denormal */
j = bbe + (P-Emin);
else
j = P + 1 - bbbits;
#endif
bb2 += j;
bd2 += j;
i = bb2 < bd2 ? bb2 : bd2;
if (i > bs2)
i = bs2;
if (i > 0) {
bb2 -= i;
bd2 -= i;
bs2 -= i;
}
if (bb5 > 0) {
bs = pow5mult(bs, bb5);
bb1 = mult(bs, bb);
Bfree(bb);
bb = bb1;
}
if (bb2 > 0)
bb = lshift(bb, bb2);
if (bd5 > 0)
bd = pow5mult(bd, bd5);
if (bd2 > 0)
bd = lshift(bd, bd2);
if (bs2 > 0)
bs = lshift(bs, bs2);
delta = diff(bb, bd);
dsign = delta->sign;
delta->sign = 0;
i = cmp(delta, bs);
if (i < 0) {
/* 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)
break;
delta = lshift(delta,Log2P);
if (cmp(delta, bs) > 0)
goto drop_down;
break;
}
if (i == 0) {
/* exactly half-way between */
if (dsign) {
if ((word0(rv) & Bndry_mask1) == Bndry_mask1
&& word1(rv) == 0xffffffff) {
/*boundary case -- increment exponent*/
word0(rv) = (word0(rv) & Exp_mask)
+ Exp_msk1
#ifdef IBM
| Exp_msk1 >> 4
#endif
;
word1(rv) = 0;
break;
}
}
else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
drop_down:
/* boundary case -- decrement exponent */
#ifdef Sudden_Underflow
L = word0(rv) & Exp_mask;
#ifdef IBM
if (L < Exp_msk1)
#else
if (L <= Exp_msk1)
#endif
goto undfl;
L -= Exp_msk1;
#else
L = (word0(rv) & Exp_mask) - Exp_msk1;
#endif
word0(rv) = L | Bndry_mask1;
word1(rv) = 0xffffffff;
#ifdef IBM
goto cont;
#else
break;
#endif
}
#ifndef ROUND_BIASED
if (!(word1(rv) & LSB))
break;
#endif
if (dsign)
rv.d += ulp(rv.d);
#ifndef ROUND_BIASED
else {
rv.d -= ulp(rv.d);
#ifndef Sudden_Underflow
if (rv.d == 0)
goto undfl;
#endif
}
#endif
break;
}
if ((aadj = ratio(delta, bs)) <= 2.) {
if (dsign)
aadj = aadj1 = 1.;
else if (word1(rv) || word0(rv) & Bndry_mask) {
#ifndef Sudden_Underflow
if (word1(rv) == Tiny1 && !word0(rv))
goto undfl;
#endif
aadj = 1.;
aadj1 = -1.;
}
else {
/* special case -- power of FLT_RADIX to be */
/* rounded down... */
if (aadj < 2./FLT_RADIX)
aadj = 1./FLT_RADIX;
else
aadj *= 0.5;
aadj1 = -aadj;
}
}
else {
aadj *= 0.5;
aadj1 = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
switch(FLT_ROUNDS) {
case 2: /* towards +infinity */
aadj1 -= 0.5;
break;
case 0: /* towards 0 */
case 3: /* towards -infinity */
aadj1 += 0.5;
}
#else
if (FLT_ROUNDS == 0)
aadj1 += 0.5;
#endif
}
y = word0(rv) & Exp_mask;
/* Check for overflow */
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
rv0.d = rv.d;
word0(rv) -= P*Exp_msk1;
adj = aadj1 * ulp(rv.d);
rv.d += adj;
if ((word0(rv) & Exp_mask) >=
Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
if (word0(rv0) == Big0 && word1(rv0) == Big1)
goto ovfl;
word0(rv) = Big0;
word1(rv) = Big1;
goto cont;
}
else
word0(rv) += P*Exp_msk1;
}
else {
#ifdef Sudden_Underflow
if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
rv0.d = rv.d;
word0(rv) += P*Exp_msk1;
adj = aadj1 * ulp(rv.d);
rv.d += adj;
#ifdef IBM
if ((word0(rv) & Exp_mask) < P*Exp_msk1)
#else
if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
#endif
{
if (word0(rv0) == Tiny0
&& word1(rv0) == Tiny1)
goto undfl;
word0(rv) = Tiny0;
word1(rv) = Tiny1;
goto cont;
}
else
word0(rv) -= P*Exp_msk1;
}
else {
adj = aadj1 * ulp(rv.d);
rv.d += adj;
}
#else
/* 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.,
* 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);
if (!dsign)
aadj1 = -aadj1;
}
adj = aadj1 * ulp(rv.d);
rv.d += adj;
#endif
}
z = word0(rv) & Exp_mask;
if (y == z) {
/* Can we stop now? */
L = aadj;
aadj -= L;
/* The tolerances below are conservative. */
if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
if (aadj < .4999999 || aadj > .5000001)
break;
}
else if (aadj < .4999999/FLT_RADIX)
break;
}
cont:
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(delta);
}
retfree:
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(bd0);
Bfree(delta);
ret:
if (se)
*se = (char *)s;
return sign ? -rv.d : rv.d;
}