Skip to content

C
cpython
Commit add28234 authored by Mark Dickinson's avatar Mark Dickinson
Browse files
Options

Merged revisions 77614-77616,77663 via svnmerge from 

svn+ssh://pythondev@svn.python.org/python/trunk

........
  r77614 | mark.dickinson | 2010-01-20 17:36:31 +0000 (Wed, 20 Jan 2010) | 5 lines

  Various dtoa.c cleanups.  1. Despagghetify _Py_dg_strtod parsing code
  and exit points.  2. Simplify bigcomp comparison loop.  3. Don't set
  ERANGE on _Py_dg_strtod underflow (it was set inconsistently anyway).
  4. Remove unused dsign field from BCinfo struct.
........
  r77615 | mark.dickinson | 2010-01-20 18:02:41 +0000 (Wed, 20 Jan 2010) | 1 line

  Don't try to put a value into a NULL pointer.
........
  r77616 | mark.dickinson | 2010-01-20 21:23:25 +0000 (Wed, 20 Jan 2010) | 1 line

  Additional explanatory comments for _Py_dg_strtod.
........
  r77663 | mark.dickinson | 2010-01-21 17:02:53 +0000 (Thu, 21 Jan 2010) | 1 line

  Additional testcases for strtod.
........
parent 577473fe
Hide whitespace changes
Inline Side-by-side
Showing with 334 additions and 189 deletions
Lib/test/test_strtod.py
View file @ add28234
...@@ -100,6 +100,49 @@ class StrtodTests(unittest.TestCase): ...@@ -100,6 +100,49 @@ class StrtodTests(unittest.TestCase):
"Incorrectly rounded str->float conversion for {}: " "Incorrectly rounded str->float conversion for {}: "
"expected {}, got {}".format(s, expected, got)) "expected {}, got {}".format(s, expected, got))
def test_short_halfway_cases(self):
# exact halfway cases with a small number of significant digits
for k in 0, 5, 10, 15, 20:
# upper = smallest integer >= 2**54/5**k
upper = -(-2**54//5**k)
# lower = smallest odd number >= 2**53/5**k
lower = -(-2**53//5**k)
if lower % 2 == 0:
lower += 1
for i in range(10 * TEST_SIZE):
# Select a random odd n in [2**53/5**k,
# 2**54/5**k). Then n * 10**k gives a halfway case
# with small number of significant digits.
n, e = random.randrange(lower, upper, 2), k
# Remove any additional powers of 5.
while n % 5 == 0:
n, e = n // 5, e + 1
assert n % 10 in (1, 3, 7, 9)
# Try numbers of the form n * 2**p2 * 10**e, p2 >= 0,
# until n * 2**p2 has more than 20 significant digits.
digits, exponent = n, e
while digits < 10**20:
s = '{}e{}'.format(digits, exponent)
self.check_strtod(s)
# Same again, but with extra trailing zeros.
s = '{}e{}'.format(digits * 10**40, exponent - 40)
self.check_strtod(s)
digits *= 2
# Try numbers of the form n * 5**p2 * 10**(e - p5), p5
# >= 0, with n * 5**p5 < 10**20.
digits, exponent = n, e
while digits < 10**20:
s = '{}e{}'.format(digits, exponent)
self.check_strtod(s)
# Same again, but with extra trailing zeros.
s = '{}e{}'.format(digits * 10**40, exponent - 40)
self.check_strtod(s)
digits *= 5
exponent -= 1
def test_halfway_cases(self): def test_halfway_cases(self):
# test halfway cases for the round-half-to-even rule # test halfway cases for the round-half-to-even rule
for i in range(1000): for i in range(1000):
...@@ -164,10 +207,10 @@ class StrtodTests(unittest.TestCase): ...@@ -164,10 +207,10 @@ class StrtodTests(unittest.TestCase):
self.check_strtod(s) self.check_strtod(s)
def test_bigcomp(self): def test_bigcomp(self):
DIG10 = 10**50 for ndigs in 5, 10, 14, 15, 16, 17, 18, 19, 20, 40, 41, 50:
for i in range(1000): dig10 = 10**ndigs
for j in range(TEST_SIZE): for i in range(100 * TEST_SIZE):
digits = random.randrange(DIG10) digits = random.randrange(dig10)
exponent = random.randrange(-400, 400) exponent = random.randrange(-400, 400)
s = '{}e{}'.format(digits, exponent) s = '{}e{}'.format(digits, exponent)
self.check_strtod(s) self.check_strtod(s)
...@@ -254,11 +297,59 @@ class StrtodTests(unittest.TestCase): ...@@ -254,11 +297,59 @@ class StrtodTests(unittest.TestCase):
# demonstration that original fix for issue 7632 bug 1 was # demonstration that original fix for issue 7632 bug 1 was
# buggy; the exit condition was too strong # buggy; the exit condition was too strong
'247032822920623295e-341', '247032822920623295e-341',
# demonstrate similar problem to issue 7632 bug1: crash
# with 'oversized quotient in quorem' message.
'99037485700245683102805043437346965248029601286431e-373',
'99617639833743863161109961162881027406769510558457e-373',
'98852915025769345295749278351563179840130565591462e-372',
'99059944827693569659153042769690930905148015876788e-373',
'98914979205069368270421829889078356254059760327101e-372',
# issue 7632 bug 5: the following 2 strings convert differently # issue 7632 bug 5: the following 2 strings convert differently
'1000000000000000000000000000000000000000e-16', '1000000000000000000000000000000000000000e-16',
'10000000000000000000000000000000000000000e-17', '10000000000000000000000000000000000000000e-17',
# issue 7632 bug 7
'991633793189150720000000000000000000000000000000000000000e-33',
# And another, similar, failing halfway case
'4106250198039490000000000000000000000000000000000000000e-38',
# issue 7632 bug 8: the following produced 10.0 # issue 7632 bug 8: the following produced 10.0
'10.900000000000000012345678912345678912345', '10.900000000000000012345678912345678912345',
# exercise exit conditions in bigcomp comparison loop
'2602129298404963083833853479113577253105939995688e2',
'260212929840496308383385347911357725310593999568896e0',
'26021292984049630838338534791135772531059399956889601e-2',
'260212929840496308383385347911357725310593999568895e0',
'260212929840496308383385347911357725310593999568897e0',
'260212929840496308383385347911357725310593999568996e0',
'260212929840496308383385347911357725310593999568866e0',
# 2**53
'9007199254740992.00',
# 2**1024 - 2**970: exact overflow boundary. All values
# smaller than this should round to something finite; any value
# greater than or equal to this one overflows.
'179769313486231580793728971405303415079934132710037' #...
'826936173778980444968292764750946649017977587207096' #...
'330286416692887910946555547851940402630657488671505' #...
'820681908902000708383676273854845817711531764475730' #...
'270069855571366959622842914819860834936475292719074' #...
'168444365510704342711559699508093042880177904174497792',
# 2**1024 - 2**970 - tiny
'179769313486231580793728971405303415079934132710037' #...
'826936173778980444968292764750946649017977587207096' #...
'330286416692887910946555547851940402630657488671505' #...
'820681908902000708383676273854845817711531764475730' #...
'270069855571366959622842914819860834936475292719074' #...
'168444365510704342711559699508093042880177904174497791.999',
# 2**1024 - 2**970 + tiny
'179769313486231580793728971405303415079934132710037' #...
'826936173778980444968292764750946649017977587207096' #...
'330286416692887910946555547851940402630657488671505' #...
'820681908902000708383676273854845817711531764475730' #...
'270069855571366959622842914819860834936475292719074' #...
'168444365510704342711559699508093042880177904174497792.001',
# 1 - 2**-54, +-tiny
'999999999999999944488848768742172978818416595458984375e-54',
'9999999999999999444888487687421729788184165954589843749999999e-54',
'9999999999999999444888487687421729788184165954589843750000001e-54',
] ]
for s in test_strings: for s in test_strings:
self.check_strtod(s) self.check_strtod(s)
......
Python/dtoa.c
View file @ add28234
...@@ -270,7 +270,7 @@ typedef union { double d; ULong L[2]; } U; ...@@ -270,7 +270,7 @@ typedef union { double d; ULong L[2]; } U;
typedef struct BCinfo BCinfo; typedef struct BCinfo BCinfo;
struct struct
BCinfo { BCinfo {
int dsign, e0, nd, nd0, scale; int e0, nd, nd0, scale;
}; };
#define FFFFFFFF 0xffffffffUL #define FFFFFFFF 0xffffffffUL
...@@ -967,8 +967,8 @@ diff(Bigint *a, Bigint *b) ...@@ -967,8 +967,8 @@ diff(Bigint *a, Bigint *b)
return c; return c;
} }
/* Given a positive normal double x, return the difference between x and the next /* Given a positive normal double x, return the difference between x and the
double up. Doesn't give correct results for subnormals. */ next double up. Doesn't give correct results for subnormals. */
static double static double
ulp(U *x) ulp(U *x)
...@@ -1276,9 +1276,6 @@ sulp(U *x, BCinfo *bc) ...@@ -1276,9 +1276,6 @@ sulp(U *x, BCinfo *bc)
bc is a struct containing information gathered during the parsing and bc is a struct containing information gathered during the parsing and
estimation steps of _Py_dg_strtod. Description of fields follows: estimation steps of _Py_dg_strtod. Description of fields follows:
bc->dsign is 1 if rv < decimal value, 0 if rv >= decimal value. In
normal use, it should almost always be 1 when bigcomp is entered.
bc->e0 gives the exponent of the input value, such that dv = (integer bc->e0 gives the exponent of the input value, such that dv = (integer
given by the bd->nd digits of s0) * 10**e0 given by the bd->nd digits of s0) * 10**e0
...@@ -1387,47 +1384,37 @@ bigcomp(U *rv, const char *s0, BCinfo *bc) ...@@ -1387,47 +1384,37 @@ bigcomp(U *rv, const char *s0, BCinfo *bc)
} }
} }
/* if b >= d, round down */ /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 ==
if (cmp(b, d) >= 0) { * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing
* a number in the range [0.1, 1). */
if (cmp(b, d) >= 0)
/* b/d >= 1 */
dd = -1; dd = -1;
goto ret; else {
} i = 0;
for(;;) {
b = multadd(b, 10, 0);
if (b == NULL) {
Bfree(d);
return -1;
}
dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
i++;
/* Compare b/d with s0 */ if (dd)
for(i = 0; i < nd0; i++) { break;
b = multadd(b, 10, 0); if (!b->x[0] && b->wds == 1) {
if (b == NULL) { /* b/d == 0 */
Bfree(d); dd = i < nd;
return -1; break;
} }
dd = *s0++ - '0' - quorem(b, d); if (!(i < nd)) {
if (dd) /* b/d != 0, but digits of s0 exhausted */
goto ret; dd = -1;
if (!b->x[0] && b->wds == 1) { break;
if (i < nd - 1) }
dd = 1;
goto ret;
}
}
s0++;
for(; i < nd; i++) {
b = multadd(b, 10, 0);
if (b == NULL) {
Bfree(d);
return -1;
}
dd = *s0++ - '0' - quorem(b, d);
if (dd)
goto ret;
if (!b->x[0] && b->wds == 1) {
if (i < nd - 1)
dd = 1;
goto ret;
} }
} }
if (b->x[0] || b->wds > 1)
dd = -1;
ret:
Bfree(b); Bfree(b);
Bfree(d); Bfree(d);
if (dd > 0 || (dd == 0 && odd)) if (dd > 0 || (dd == 0 && odd))
...@@ -1438,128 +1425,130 @@ bigcomp(U *rv, const char *s0, BCinfo *bc) ...@@ -1438,128 +1425,130 @@ bigcomp(U *rv, const char *s0, BCinfo *bc)
double double
_Py_dg_strtod(const char *s00, char **se) _Py_dg_strtod(const char *s00, char **se)
{ {
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1, error; int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e, e1, error;
int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; int esign, i, j, k, lz, nd, nd0, sign;
const char *s, *s0, *s1; const char *s, *s0, *s1;
double aadj, aadj1; double aadj, aadj1;
U aadj2, adj, rv, rv0; U aadj2, adj, rv, rv0;
ULong y, z, abse; ULong y, z, abs_exp;
Long L; Long L;
BCinfo bc; BCinfo bc;
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
sign = nz0 = nz = 0;
dval(&rv) = 0.; dval(&rv) = 0.;
for(s = s00;;s++) switch(*s) {
case '-': /* Start parsing. */
sign = 1; c = *(s = s00);
/* no break */
case '+': /* Parse optional sign, if present. */
if (*++s) sign = 0;
goto break2; switch (c) {
/* no break */ case '-':
case 0: sign = 1;
goto ret0; /* no break */
/* modify original dtoa.c so that it doesn't accept leading whitespace case '+':
case '\t': c = *++s;
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;
for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) /* Skip leading zeros: lz is true iff there were leading zeros. */
; s1 = s;
nd0 = nd; while (c == '0')
c = *++s;
lz = s != s1;
/* Point s0 at the first nonzero digit (if any). nd0 will be the position
of the point relative to s0. nd will be the total number of digits
ignoring leading zeros. */
s0 = s1 = s;
while ('0' <= c && c <= '9')
c = *++s;
nd0 = nd = s - s1;
/* Parse decimal point and following digits. */
if (c == '.') { if (c == '.') {
c = *++s; c = *++s;
if (!nd) { if (!nd) {
for(; c == '0'; c = *++s) s1 = s;
nz++; while (c == '0')
if (c > '0' && c <= '9') { c = *++s;
s0 = s; lz = lz || s != s1;
nf += nz; nd0 -= s - s1;
nz = 0; s0 = s;
goto have_dig;
}
goto dig_done;
}
for(; c >= '0' && c <= '9'; c = *++s) {
have_dig:
nz++;
if (c -= '0') {
nf += nz;
nd += nz;
nz = 0;
}
} }
s1 = s;
while ('0' <= c && c <= '9')
c = *++s;
nd += s - s1;
}
/* Now lz is true if and only if there were leading zero digits, and nd
gives the total number of digits ignoring leading zeros. A valid input
must have at least one digit. */
if (!nd && !lz) {
if (se)
*se = (char *)s00;
goto parse_error;
} }
dig_done:
/* Parse exponent. */
e = 0; e = 0;
if (c == 'e' || c == 'E') { if (c == 'e' || c == 'E') {
if (!nd && !nz && !nz0) {
goto ret0;
}
s00 = s; s00 = s;
c = *++s;
/* Exponent sign. */
esign = 0; esign = 0;
switch(c = *++s) { switch (c) {
case '-': case '-':
esign = 1; esign = 1;
/* no break */
case '+': case '+':
c = *++s; c = *++s;
} }
if (c >= '0' && c <= '9') {
while(c == '0') /* Skip zeros. lz is true iff there are leading zeros. */
c = *++s; s1 = s;
if (c > '0' && c <= '9') { while (c == '0')
abse = c - '0'; c = *++s;
s1 = s; lz = s != s1;
while((c = *++s) >= '0' && c <= '9')
abse = 10*abse + c - '0'; /* Get absolute value of the exponent. */
if (s - s1 > 8 || abse > MAX_ABS_EXP) s1 = s;
/* Avoid confusion from exponents abs_exp = 0;
* so large that e might overflow. while ('0' <= c && c <= '9') {
*/ abs_exp = 10*abs_exp + (c - '0');
e = (int)MAX_ABS_EXP; /* safe for 16 bit ints */ c = *++s;
else
e = (int)abse;
if (esign)
e = -e;
}
else
e = 0;
} }
/* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if
there are at most 9 significant exponent digits then overflow is
impossible. */
if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
e = (int)MAX_ABS_EXP;
else else
e = (int)abs_exp;
if (esign)
e = -e;
/* A valid exponent must have at least one digit. */
if (s == s1 && !lz)
s = s00; s = s00;
} }
if (!nd) {
if (!nz && !nz0) { /* Adjust exponent to take into account position of the point. */
ret0: e -= nd - nd0;
s = s00; if (nd0 <= 0)
sign = 0;
}
goto ret;
}
e -= nf;
if (!nd0)
nd0 = nd; nd0 = nd;
/* strip trailing zeros */ /* Finished parsing. Set se to indicate how far we parsed */
if (se)
*se = (char *)s;
/* If all digits were zero, exit with return value +-0.0. Otherwise,
strip trailing zeros: scan back until we hit a nonzero digit. */
if (!nd)
goto ret;
for (i = nd; i > 0; ) { for (i = nd; i > 0; ) {
/* scan back until we hit a nonzero digit. significant digit 'i'
is s0[i] if i < nd0, s0[i+1] if i >= nd0. */
--i; --i;
if (s0[i < nd0 ? i : i+1] != '0') { if (s0[i < nd0 ? i : i+1] != '0') {
++i; ++i;
...@@ -1571,28 +1560,21 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1571,28 +1560,21 @@ _Py_dg_strtod(const char *s00, char **se)
if (nd0 > nd) if (nd0 > nd)
nd0 = nd; nd0 = nd;
/* Now we have nd0 digits, starting at s0, followed by a /* Summary of parsing results. After parsing, and dealing with zero
* decimal point, followed by nd-nd0 digits. The number we're * inputs, we have values s0, nd0, nd, e, sign, where:
* after is the integer represented by those digits times
* 10**e */
bc.e0 = e1 = e;
/* Summary of parsing results. The parsing stage gives values
* s0, nd0, nd, e, sign, where:
* *
* - s0 points to the first significant digit of the input string s00; * - s0 points to the first significant digit of the input string
* *
* - nd is the total number of significant digits (here, and * - nd is the total number of significant digits (here, and
* below, 'significant digits' means the set of digits of the * below, 'significant digits' means the set of digits of the
* significand of the input that remain after ignoring leading * significand of the input that remain after ignoring leading
* and trailing zeros. * and trailing zeros).
* *
* - nd0 indicates the position of the decimal point (if * - nd0 indicates the position of the decimal point, if present; it
* present): so the nd significant digits are in s0[0:nd0] and * satisfies 1 <= nd0 <= nd. The nd significant digits are in
* s0[nd0+1:nd+1] using the usual Python half-open slice * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice
* notation. (If nd0 < nd, then s0[nd0] necessarily contains * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if
* a '.' character; if nd0 == nd, then it could be anything.) * nd0 == nd, then s0[nd0] could be any non-digit character.)
* *
* - e is the adjusted exponent: the absolute value of the number * - e is the adjusted exponent: the absolute value of the number
* represented by the original input string is n * 10**e, where * represented by the original input string is n * 10**e, where
...@@ -1614,6 +1596,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1614,6 +1596,7 @@ _Py_dg_strtod(const char *s00, char **se)
* gives the value represented by the first min(16, nd) sig. digits. * gives the value represented by the first min(16, nd) sig. digits.
*/ */
bc.e0 = e1 = e;
y = z = 0; y = z = 0;
for (i = 0; i < nd; i++) { for (i = 0; i < nd; i++) {
if (i < 9) if (i < 9)
...@@ -1666,14 +1649,8 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1666,14 +1649,8 @@ _Py_dg_strtod(const char *s00, char **se)
if ((i = e1 & 15)) if ((i = e1 & 15))
dval(&rv) *= tens[i]; dval(&rv) *= tens[i];
if (e1 &= ~15) { if (e1 &= ~15) {
if (e1 > DBL_MAX_10_EXP) { if (e1 > DBL_MAX_10_EXP)
ovfl: goto ovfl;
errno = ERANGE;
/* Can't trust HUGE_VAL */
word0(&rv) = Exp_mask;
word1(&rv) = 0;
goto ret;
}
e1 >>= 4; e1 >>= 4;
for(j = 0; e1 > 1; j++, e1 >>= 1) for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1) if (e1 & 1)
...@@ -1695,6 +1672,16 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1695,6 +1672,16 @@ _Py_dg_strtod(const char *s00, char **se)
} }
} }
else if (e1 < 0) { else if (e1 < 0) {
/* The input decimal value lies in [10**e1, 10**(e1+16)).
If e1 <= -512, underflow immediately.
If e1 <= -256, set bc.scale to 2*P.
So for input value < 1e-256, bc.scale is always set;
for input value >= 1e-240, bc.scale is never set.
For input values in [1e-256, 1e-240), bc.scale may or may
not be set. */
e1 = -e1; e1 = -e1;
if ((i = e1 & 15)) if ((i = e1 & 15))
dval(&rv) /= tens[i]; dval(&rv) /= tens[i];
...@@ -1719,12 +1706,8 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1719,12 +1706,8 @@ _Py_dg_strtod(const char *s00, char **se)
else else
word1(&rv) &= 0xffffffff << j; word1(&rv) &= 0xffffffff << j;
} }
if (!dval(&rv)) { if (!dval(&rv))
undfl: goto undfl;
dval(&rv) = 0.;
errno = ERANGE;
goto ret;
}
} }
} }
...@@ -1769,7 +1752,34 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1769,7 +1752,34 @@ _Py_dg_strtod(const char *s00, char **se)
if (bd0 == NULL) if (bd0 == NULL)
goto failed_malloc; goto failed_malloc;
/* Notation for the comments below. Write:
- dv for the absolute value of the number represented by the original
decimal input string.
- if we've truncated dv, write tdv for the truncated value.
Otherwise, set tdv == dv.
- srv for the quantity rv/2^bc.scale; so srv is the current binary
approximation to tdv (and dv). It should be exactly representable
in an IEEE 754 double.
*/
for(;;) { for(;;) {
/* This is the main correction loop for _Py_dg_strtod.
We've got a decimal value tdv, and a floating-point approximation
srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is
close enough (i.e., within 0.5 ulps) to tdv, and to compute a new
approximation if not.
To determine whether srv is close enough to tdv, compute integers
bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv)
respectively, and then use integer arithmetic to determine whether
|tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv).
*/
bd = Balloc(bd0->k); bd = Balloc(bd0->k);
if (bd == NULL) { if (bd == NULL) {
Bfree(bd0); Bfree(bd0);
...@@ -1782,6 +1792,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1782,6 +1792,7 @@ _Py_dg_strtod(const char *s00, char **se)
Bfree(bd0); Bfree(bd0);
goto failed_malloc; goto failed_malloc;
} }
/* tdv = bd * 10^e; srv = bb * 2^(bbe - scale) */
bs = i2b(1); bs = i2b(1);
if (bs == NULL) { if (bs == NULL) {
Bfree(bb); Bfree(bb);
...@@ -1802,6 +1813,17 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1802,6 +1813,17 @@ _Py_dg_strtod(const char *s00, char **se)
bb2 += bbe; bb2 += bbe;
else else
bd2 -= bbe; bd2 -= bbe;
/* At this stage e = bd2 - bb2 + bbe = bd5 - bb5, so:
tdv = bd * 2^(bbe + bd2 - bb2) * 5^(bd5 - bb5)
srv = bb * 2^(bbe - scale).
Compute j such that
0.5 ulp(srv) = 2^(bbe - scale - j)
*/
bs2 = bb2; bs2 = bb2;
j = bbe - bc.scale; j = bbe - bc.scale;
i = j + bbbits - 1; /* logb(rv) */ i = j + bbbits - 1; /* logb(rv) */
...@@ -1809,9 +1831,26 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1809,9 +1831,26 @@ _Py_dg_strtod(const char *s00, char **se)
j += P - Emin; j += P - Emin;
else else
j = P + 1 - bbbits; j = P + 1 - bbbits;
/* Now we have:
M * tdv = bd * 2^(bd2 + scale + j) * 5^bd5
M * srv = bb * 2^(bb2 + j) * 5^bb5
M * 0.5 ulp(srv) = 2^bs2 * 5^bb5
where M is the constant (currently) represented by 2^(j + scale +
bb2 - bbe) * 5^bb5.
*/
bb2 += j; bb2 += j;
bd2 += j; bd2 += j;
bd2 += bc.scale; bd2 += bc.scale;
/* After the adjustments above, tdv, srv and 0.5 ulp(srv) are
proportional to: bd * 2^bd2 * 5^bd5, bb * 2^bb2 * 5^bb5, and
bs * 2^bs2 * 5^bb5, respectively. */
/* Remove excess powers of 2. i = min(bb2, bd2, bs2). */
i = bb2 < bd2 ? bb2 : bd2; i = bb2 < bd2 ? bb2 : bd2;
if (i > bs2) if (i > bs2)
i = bs2; i = bs2;
...@@ -1820,6 +1859,8 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1820,6 +1859,8 @@ _Py_dg_strtod(const char *s00, char **se)
bd2 -= i; bd2 -= i;
bs2 -= i; bs2 -= i;
} }
/* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
if (bb5 > 0) { if (bb5 > 0) {
bs = pow5mult(bs, bb5); bs = pow5mult(bs, bb5);
if (bs == NULL) { if (bs == NULL) {
...@@ -1874,6 +1915,11 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1874,6 +1915,11 @@ _Py_dg_strtod(const char *s00, char **se)
goto failed_malloc; goto failed_malloc;
} }
} }
/* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv),
respectively. Compute the difference |tdv - srv|, and compare
with 0.5 ulp(srv). */
delta = diff(bb, bd); delta = diff(bb, bd);
if (delta == NULL) { if (delta == NULL) {
Bfree(bb); Bfree(bb);
...@@ -1882,11 +1928,11 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1882,11 +1928,11 @@ _Py_dg_strtod(const char *s00, char **se)
Bfree(bd0); Bfree(bd0);
goto failed_malloc; goto failed_malloc;
} }
bc.dsign = delta->sign; dsign = delta->sign;
delta->sign = 0; delta->sign = 0;
i = cmp(delta, bs); i = cmp(delta, bs);
if (bc.nd > nd && i <= 0) { if (bc.nd > nd && i <= 0) {
if (bc.dsign) if (dsign)
break; /* Must use bigcomp(). */ break; /* Must use bigcomp(). */
/* Here rv overestimates the truncated decimal value by at most /* Here rv overestimates the truncated decimal value by at most
...@@ -1908,7 +1954,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1908,7 +1954,7 @@ _Py_dg_strtod(const char *s00, char **se)
rv / 2^bc.scale >= 2^-1021. */ rv / 2^bc.scale >= 2^-1021. */
if (j - bc.scale >= 2) { if (j - bc.scale >= 2) {
dval(&rv) -= 0.5 * sulp(&rv, &bc); dval(&rv) -= 0.5 * sulp(&rv, &bc);
break; break; /* Use bigcomp. */
} }
} }
...@@ -1922,7 +1968,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1922,7 +1968,7 @@ _Py_dg_strtod(const char *s00, char **se)
/* Error is less than half an ulp -- check for /* Error is less than half an ulp -- check for
* special case of mantissa a power of two. * special case of mantissa a power of two.
*/ */
if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
|| (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
) { ) {
break; break;
...@@ -1945,7 +1991,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1945,7 +1991,7 @@ _Py_dg_strtod(const char *s00, char **se)
} }
if (i == 0) { if (i == 0) {
/* exactly half-way between */ /* exactly half-way between */
if (bc.dsign) { if (dsign) {
if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
&& word1(&rv) == ( && word1(&rv) == (
(bc.scale && (bc.scale &&
...@@ -1957,7 +2003,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1957,7 +2003,7 @@ _Py_dg_strtod(const char *s00, char **se)
+ Exp_msk1 + Exp_msk1
; ;
word1(&rv) = 0; word1(&rv) = 0;
bc.dsign = 0; dsign = 0;
break; break;
} }
} }
...@@ -1972,7 +2018,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1972,7 +2018,7 @@ _Py_dg_strtod(const char *s00, char **se)
/* accept rv */ /* accept rv */
break; break;
/* rv = smallest denormal */ /* rv = smallest denormal */
if (bc.nd >nd) if (bc.nd > nd)
break; break;
goto undfl; goto undfl;
} }
...@@ -1984,7 +2030,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1984,7 +2030,7 @@ _Py_dg_strtod(const char *s00, char **se)
} }
if (!(word1(&rv) & LSB)) if (!(word1(&rv) & LSB))
break; break;
if (bc.dsign) if (dsign)
dval(&rv) += ulp(&rv); dval(&rv) += ulp(&rv);
else { else {
dval(&rv) -= ulp(&rv); dval(&rv) -= ulp(&rv);
...@@ -1994,11 +2040,11 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -1994,11 +2040,11 @@ _Py_dg_strtod(const char *s00, char **se)
goto undfl; goto undfl;
} }
} }
bc.dsign = 1 - bc.dsign; dsign = 1 - dsign;
break; break;
} }
if ((aadj = ratio(delta, bs)) <= 2.) { if ((aadj = ratio(delta, bs)) <= 2.) {
if (bc.dsign) if (dsign)
aadj = aadj1 = 1.; aadj = aadj1 = 1.;
else if (word1(&rv) || word0(&rv) & Bndry_mask) { else if (word1(&rv) || word0(&rv) & Bndry_mask) {
if (word1(&rv) == Tiny1 && !word0(&rv)) { if (word1(&rv) == Tiny1 && !word0(&rv)) {
...@@ -2022,7 +2068,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -2022,7 +2068,7 @@ _Py_dg_strtod(const char *s00, char **se)
} }
else { else {
aadj *= 0.5; aadj *= 0.5;
aadj1 = bc.dsign ? aadj : -aadj; aadj1 = dsign ? aadj : -aadj;
if (Flt_Rounds == 0) if (Flt_Rounds == 0)
aadj1 += 0.5; aadj1 += 0.5;
} }
...@@ -2058,7 +2104,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -2058,7 +2104,7 @@ _Py_dg_strtod(const char *s00, char **se)
if ((z = (ULong)aadj) <= 0) if ((z = (ULong)aadj) <= 0)
z = 1; z = 1;
aadj = z; aadj = z;
aadj1 = bc.dsign ? aadj : -aadj; aadj1 = dsign ? aadj : -aadj;
} }
dval(&aadj2) = aadj1; dval(&aadj2) = aadj1;
word0(&aadj2) += (2*P+1)*Exp_msk1 - y; word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
...@@ -2075,7 +2121,7 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -2075,7 +2121,7 @@ _Py_dg_strtod(const char *s00, char **se)
L = (Long)aadj; L = (Long)aadj;
aadj -= L; aadj -= L;
/* The tolerances below are conservative. */ /* The tolerances below are conservative. */
if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
if (aadj < .4999999 || aadj > .5000001) if (aadj < .4999999 || aadj > .5000001)
break; break;
} }
...@@ -2104,20 +2150,28 @@ _Py_dg_strtod(const char *s00, char **se) ...@@ -2104,20 +2150,28 @@ _Py_dg_strtod(const char *s00, char **se)
word0(&rv0) = Exp_1 - 2*P*Exp_msk1; word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
word1(&rv0) = 0; word1(&rv0) = 0;
dval(&rv) *= dval(&rv0); dval(&rv) *= dval(&rv0);
/* try to avoid the bug of testing an 8087 register value */
if (!(word0(&rv) & Exp_mask))
errno = ERANGE;
} }
ret: ret:
if (se)
*se = (char *)s;
return sign ? -dval(&rv) : dval(&rv); return sign ? -dval(&rv) : dval(&rv);
parse_error:
return 0.0;
failed_malloc: failed_malloc:
if (se)
*se = (char *)s00;
errno = ENOMEM; errno = ENOMEM;
return -1.0; return -1.0;
undfl:
return sign ? -0.0 : 0.0;
ovfl:
errno = ERANGE;
/* Can't trust HUGE_VAL */
word0(&rv) = Exp_mask;
word1(&rv) = 0;
return sign ? -dval(&rv) : dval(&rv);
} }
static char * static char *
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment
GitLab Nexedi Edition | About GitLab | About Nexedi | 沪ICP备2021021310号-2 | 沪ICP备2021021310号-7