From 9000c1614d183345165271c42926da80102a2a36 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Thu, 14 Jan 2010 15:43:57 +0000 Subject: Merged revisions 77494 via svnmerge from svn+ssh://pythondev@svn.python.org/python/branches/py3k ................ r77494 | mark.dickinson | 2010-01-14 15:37:49 +0000 (Thu, 14 Jan 2010) | 41 lines Merged revisions 77477-77478,77481-77483,77490-77493 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r77477 | mark.dickinson | 2010-01-13 18:21:53 +0000 (Wed, 13 Jan 2010) | 1 line Add comments explaining the role of the bigcomp function in dtoa.c. ........ r77478 | mark.dickinson | 2010-01-13 19:02:37 +0000 (Wed, 13 Jan 2010) | 1 line Clarify that sulp expects a nonnegative input, but that +0.0 is fine. ........ r77481 | mark.dickinson | 2010-01-13 20:55:03 +0000 (Wed, 13 Jan 2010) | 1 line Simplify and annotate the bigcomp function, removing unused special cases. ........ r77482 | mark.dickinson | 2010-01-13 22:15:53 +0000 (Wed, 13 Jan 2010) | 1 line Fix buggy comparison: LHS of comparison was being treated as unsigned. ........ r77483 | mark.dickinson | 2010-01-13 22:20:10 +0000 (Wed, 13 Jan 2010) | 1 line More dtoa.c cleanup; remove the need for bc.dplen, bc.dp0 and bc.dp1. ........ r77490 | mark.dickinson | 2010-01-14 13:02:36 +0000 (Thu, 14 Jan 2010) | 1 line Fix off-by-one error introduced in r77483. I have a test for this, but it currently fails due to a different dtoa.c bug; I'll add the test once that bug is fixed. ........ r77491 | mark.dickinson | 2010-01-14 13:14:49 +0000 (Thu, 14 Jan 2010) | 1 line Issue 7632: fix a dtoa.c bug (bug 6) causing incorrect rounding. Tests to follow. ........ r77492 | mark.dickinson | 2010-01-14 14:40:20 +0000 (Thu, 14 Jan 2010) | 1 line Issue 7632: fix incorrect rounding for long input strings with values very close to a power of 2. (See Bug 4 in the tracker discussion.) ........ r77493 | mark.dickinson | 2010-01-14 15:22:33 +0000 (Thu, 14 Jan 2010) | 1 line Issue #7632: add tests for bugs fixed so far. ........ ................ --- Lib/test/test_strtod.py | 269 ++++++++++++++++++++++++++++++++++++++++++++++++ Python/dtoa.c | 258 ++++++++++++++++++++++++++-------------------- 2 files changed, 414 insertions(+), 113 deletions(-) create mode 100644 Lib/test/test_strtod.py diff --git a/Lib/test/test_strtod.py b/Lib/test/test_strtod.py new file mode 100644 index 0000000..79cfc88 --- /dev/null +++ b/Lib/test/test_strtod.py @@ -0,0 +1,269 @@ +# Tests for the correctly-rounded string -> float conversions +# introduced in Python 2.7 and 3.1. + +import random +import struct +import unittest +import re +import sys +import test.support + +# Correctly rounded str -> float in pure Python, for comparison. + +strtod_parser = re.compile(r""" # A numeric string consists of: + (?P[-+])? # an optional sign, followed by + (?=\d|\.\d) # a number with at least one digit + (?P\d*) # having a (possibly empty) integer part + (?:\.(?P\d*))? # followed by an optional fractional part + (?:E(?P[-+]?\d+))? # and an optional exponent + \Z +""", re.VERBOSE | re.IGNORECASE).match + +def strtod(s, mant_dig=53, min_exp = -1021, max_exp = 1024): + """Convert a finite decimal string to a hex string representing an + IEEE 754 binary64 float. Return 'inf' or '-inf' on overflow. + This function makes no use of floating-point arithmetic at any + stage.""" + + # parse string into a pair of integers 'a' and 'b' such that + # abs(decimal value) = a/b, along with a boolean 'negative'. + m = strtod_parser(s) + if m is None: + raise ValueError('invalid numeric string') + fraction = m.group('frac') or '' + intpart = int(m.group('int') + fraction) + exp = int(m.group('exp') or '0') - len(fraction) + negative = m.group('sign') == '-' + a, b = intpart*10**max(exp, 0), 10**max(0, -exp) + + # quick return for zeros + if not a: + return '-0x0.0p+0' if negative else '0x0.0p+0' + + # compute exponent e for result; may be one too small in the case + # that the rounded value of a/b lies in a different binade from a/b + d = a.bit_length() - b.bit_length() + d += (a >> d if d >= 0 else a << -d) >= b + e = max(d, min_exp) - mant_dig + + # approximate a/b by number of the form q * 2**e; adjust e if necessary + a, b = a << max(-e, 0), b << max(e, 0) + q, r = divmod(a, b) + if 2*r > b or 2*r == b and q & 1: + q += 1 + if q.bit_length() == mant_dig+1: + q //= 2 + e += 1 + + # double check that (q, e) has the right form + assert q.bit_length() <= mant_dig and e >= min_exp - mant_dig + assert q.bit_length() == mant_dig or e == min_exp - mant_dig + + # check for overflow and underflow + if e + q.bit_length() > max_exp: + return '-inf' if negative else 'inf' + if not q: + return '-0x0.0p+0' if negative else '0x0.0p+0' + + # for hex representation, shift so # bits after point is a multiple of 4 + hexdigs = 1 + (mant_dig-2)//4 + shift = 3 - (mant_dig-2)%4 + q, e = q << shift, e - shift + return '{}0x{:x}.{:0{}x}p{:+d}'.format( + '-' if negative else '', + q // 16**hexdigs, + q % 16**hexdigs, + hexdigs, + e + 4*hexdigs) + +TEST_SIZE = 10 + +@unittest.skipUnless(getattr(sys, 'float_repr_style', '') == 'short', + "applies only when using short float repr style") +class StrtodTests(unittest.TestCase): + def check_strtod(self, s): + """Compare the result of Python's builtin correctly rounded + string->float conversion (using float) to a pure Python + correctly rounded string->float implementation. Fail if the + two methods give different results.""" + + try: + fs = float(s) + except OverflowError: + got = '-inf' if s[0] == '-' else 'inf' + else: + got = fs.hex() + expected = strtod(s) + self.assertEqual(expected, got, + "Incorrectly rounded str->float conversion for {}: " + "expected {}, got {}".format(s, expected, got)) + + def test_halfway_cases(self): + # test halfway cases for the round-half-to-even rule + for i in range(1000): + for j in range(TEST_SIZE): + # bit pattern for a random finite positive (or +0.0) float + bits = random.randrange(2047*2**52) + + # convert bit pattern to a number of the form m * 2**e + e, m = divmod(bits, 2**52) + if e: + m, e = m + 2**52, e - 1 + e -= 1074 + + # add 0.5 ulps + m, e = 2*m + 1, e - 1 + + # convert to a decimal string + if e >= 0: + digits = m << e + exponent = 0 + else: + # m * 2**e = (m * 5**-e) * 10**e + digits = m * 5**-e + exponent = e + s = '{}e{}'.format(digits, exponent) + + # for the moment, ignore errors from trailing zeros + if digits % 10 == 0: + continue + self.check_strtod(s) + + # get expected answer via struct, to triple check + #fs = struct.unpack('scale != 0 (issue 7632, bug 6). + '104308485241983990666713401708072175773165034278685' #... + '682646111762292409330928739751702404658197872319129' #... + '036519947435319418387839758990478549477777586673075' #... + '945844895981012024387992135617064532141489278815239' #... + '849108105951619997829153633535314849999674266169258' #... + '928940692239684771590065027025835804863585454872499' #... + '320500023126142553932654370362024104462255244034053' #... + '203998964360882487378334860197725139151265590832887' #... + '433736189468858614521708567646743455601905935595381' #... + '852723723645799866672558576993978025033590728687206' #... + '296379801363024094048327273913079612469982585674824' #... + '156000783167963081616214710691759864332339239688734' #... + '656548790656486646106983450809073750535624894296242' #... + '072010195710276073042036425579852459556183541199012' #... + '652571123898996574563824424330960027873516082763671875e-1075', + # demonstration that original fix for issue 7632 bug 1 was + # buggy; the exit condition was too strong + '247032822920623295e-341', + # issue 7632 bug 5: the following 2 strings convert differently + '1000000000000000000000000000000000000000e-16', + #'10000000000000000000000000000000000000000e-17', + ] + for s in test_strings: + self.check_strtod(s) + +def test_main(): + test.support.run_unittest(StrtodTests) + +if __name__ == "__main__": + test_main() diff --git a/Python/dtoa.c b/Python/dtoa.c index 1fe20f4..51895c7 100644 --- a/Python/dtoa.c +++ b/Python/dtoa.c @@ -270,7 +270,7 @@ typedef union { double d; ULong L[2]; } U; typedef struct BCinfo BCinfo; struct BCinfo { - int dp0, dp1, dplen, dsign, e0, nd, nd0, scale; + int dsign, e0, nd, nd0, scale; }; #define FFFFFFFF 0xffffffffUL @@ -437,7 +437,7 @@ multadd(Bigint *b, int m, int a) /* multiply by m and add a */ NULL on failure. */ static Bigint * -s2b(const char *s, int nd0, int nd, ULong y9, int dplen) +s2b(const char *s, int nd0, int nd, ULong y9) { Bigint *b; int i, k; @@ -451,18 +451,16 @@ s2b(const char *s, int nd0, int nd, ULong y9, int dplen) b->x[0] = y9; b->wds = 1; - i = 9; - if (9 < nd0) { - s += 9; - do { - b = multadd(b, 10, *s++ - '0'); - if (b == NULL) - return NULL; - } while(++i < nd0); - s += dplen; + if (nd <= 9) + return b; + + s += 9; + for (i = 9; i < nd0; i++) { + b = multadd(b, 10, *s++ - '0'); + if (b == NULL) + return NULL; } - else - s += dplen + 9; + s++; for(; i < nd; i++) { b = multadd(b, 10, *s++ - '0'); if (b == NULL) @@ -1130,76 +1128,120 @@ quorem(Bigint *b, Bigint *S) return q; } -/* version of ulp(x) that takes bc.scale into account. +/* sulp(x) is a version of ulp(x) that takes bc.scale into account. - Assuming that x is finite and nonzero, and x / 2^bc.scale is exactly - representable as a double, sulp(x) is equivalent to 2^bc.scale * ulp(x / - 2^bc.scale). */ + Assuming that x is finite and nonnegative (positive zero is fine + here) and x / 2^bc.scale is exactly representable as a double, + sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */ static double sulp(U *x, BCinfo *bc) { U u; - if (bc->scale && 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift) > 0) { + if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) { /* rv/2^bc->scale is subnormal */ word0(&u) = (P+2)*Exp_msk1; word1(&u) = 0; return u.d; } - else + else { + assert(word0(x) || word1(x)); /* x != 0.0 */ return ulp(x); + } } -/* return 0 on success, -1 on failure */ +/* The bigcomp function handles some hard cases for strtod, for inputs + with more than STRTOD_DIGLIM digits. It's called once an initial + estimate for the double corresponding to the input string has + already been obtained by the code in _Py_dg_strtod. + + The bigcomp function is only called after _Py_dg_strtod has found a + double value rv such that either rv or rv + 1ulp represents the + correctly rounded value corresponding to the original string. It + determines which of these two values is the correct one by + computing the decimal digits of rv + 0.5ulp and comparing them with + the corresponding digits of s0. + + In the following, write dv for the absolute value of the number represented + by the input string. + + Inputs: + + s0 points to the first significant digit of the input string. + + rv is a (possibly scaled) estimate for the closest double value to the + value represented by the original input to _Py_dg_strtod. If + bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to + the input value. + + bc is a struct containing information gathered during the parsing and + 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 + given by the bd->nd digits of s0) * 10**e0 + + bc->nd gives the total number of significant digits of s0. It will + be at least 1. + + bc->nd0 gives the number of significant digits of s0 before the + decimal separator. If there's no decimal separator, bc->nd0 == + bc->nd. + + bc->scale is the value used to scale rv to avoid doing arithmetic with + subnormal values. It's either 0 or 2*P (=106). + + Outputs: + + On successful exit, rv/2^(bc->scale) is the closest double to dv. + + Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */ static int bigcomp(U *rv, const char *s0, BCinfo *bc) { Bigint *b, *d; - int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; + int b2, bbits, d2, dd, i, nd, nd0, odd, p2, p5; - dsign = bc->dsign; + dd = 0; /* silence compiler warning about possibly unused variable */ nd = bc->nd; nd0 = bc->nd0; p5 = nd + bc->e0; - speccase = 0; - if (rv->d == 0.) { /* special case: value near underflow-to-zero */ - /* threshold was rounded to zero */ - b = i2b(1); + if (rv->d == 0.) { + /* special case because d2b doesn't handle 0.0 */ + b = i2b(0); if (b == NULL) return -1; - p2 = Emin - P + 1; - bbits = 1; - word0(rv) = (P+2) << Exp_shift; - i = 0; - { - speccase = 1; - --p2; - dsign = 0; - goto have_i; - } + p2 = Emin - P + 1; /* = -1074 for IEEE 754 binary64 */ + bbits = 0; } - else - { + else { b = d2b(rv, &p2, &bbits); if (b == NULL) return -1; + p2 -= bc->scale; } - p2 -= bc->scale; - /* floor(log2(rv)) == bbits - 1 + p2 */ - /* Check for denormal case. */ + /* now rv/2^(bc->scale) = b * 2**p2, and b has bbits significant bits */ + + /* Replace (b, p2) by (b << i, p2 - i), with i the largest integer such + that b << i has at most P significant bits and p2 - i >= Emin - P + + 1. */ i = P - bbits; - if (i > (j = P - Emin - 1 + p2)) { - i = j; - } - { - b = lshift(b, ++i); - if (b == NULL) - return -1; - b->x[0] |= 1; - } - have_i: + if (i > p2 - (Emin - P + 1)) + i = p2 - (Emin - P + 1); + /* increment i so that we shift b by an extra bit; then or-ing a 1 into + the lsb of b gives us rv/2^(bc->scale) + 0.5ulp. */ + b = lshift(b, ++i); + if (b == NULL) + return -1; + /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway + case, this is used for round to even. */ + odd = b->x[0] & 2; + b->x[0] |= 1; + p2 -= p5 + i; d = i2b(1); if (d == NULL) { @@ -1247,92 +1289,58 @@ bigcomp(U *rv, const char *s0, BCinfo *bc) } } - /* Now 10*b/d = exactly half-way between the two floating-point values - on either side of the input string. If b >= d, round down. */ + /* if b >= d, round down */ if (cmp(b, d) >= 0) { dd = -1; goto ret; } - - /* Compute first digit of 10*b/d. */ - b = multadd(b, 10, 0); - if (b == NULL) { - Bfree(d); - return -1; - } - dig = quorem(b, d); - assert(dig < 10); /* Compare b/d with s0 */ - - assert(nd > 0); - dd = 9999; /* silence gcc compiler warning */ - for(i = 0; i < nd0; ) { - if ((dd = s0[i++] - '0' - dig)) - goto ret; - if (!b->x[0] && b->wds == 1) { - if (i < nd) - dd = 1; - goto ret; - } + for(i = 0; i < nd0; i++) { b = multadd(b, 10, 0); if (b == NULL) { Bfree(d); return -1; } - dig = quorem(b,d); - } - for(j = bc->dp1; i++ < nd;) { - if ((dd = s0[j++] - '0' - dig)) + dd = *s0++ - '0' - quorem(b, d); + if (dd) goto ret; if (!b->x[0] && b->wds == 1) { - if (i < nd) + 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; } - dig = quorem(b,d); + 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(d); - if (speccase) { - if (dd <= 0) - rv->d = 0.; - } - else if (dd < 0) { - if (!dsign) /* does not happen for round-near */ - retlow1: - dval(rv) -= sulp(rv, bc); - } - else if (dd > 0) { - if (dsign) { - rethi1: - dval(rv) += sulp(rv, bc); - } - } - else { - /* Exact half-way case: apply round-even rule. */ - if (word1(rv) & 1) { - if (dsign) - goto rethi1; - goto retlow1; - } - } - + if (dd > 0 || (dd == 0 && odd)) + dval(rv) += sulp(rv, bc); return 0; } double _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, dp0, dp1, dplen, e, e1, error; int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; const char *s, *s0, *s1; double aadj, aadj1; @@ -1341,7 +1349,7 @@ _Py_dg_strtod(const char *s00, char **se) BCinfo bc; Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; - sign = nz0 = nz = bc.dplen = 0; + sign = nz0 = nz = dplen = 0; dval(&rv) = 0.; for(s = s00;;s++) switch(*s) { case '-': @@ -1380,11 +1388,11 @@ _Py_dg_strtod(const char *s00, char **se) else if (nd < 16) z = 10*z + c - '0'; nd0 = nd; - bc.dp0 = bc.dp1 = s - s0; + dp0 = dp1 = s - s0; if (c == '.') { c = *++s; - bc.dp1 = s - s0; - bc.dplen = bc.dp1 - bc.dp0; + dp1 = s - s0; + dplen = 1; if (!nd) { for(; c == '0'; c = *++s) nz++; @@ -1587,10 +1595,10 @@ _Py_dg_strtod(const char *s00, char **se) /* in IEEE arithmetic. */ i = j = 18; if (i > nd0) - j += bc.dplen; + j += dplen; for(;;) { - if (--j <= bc.dp1 && j >= bc.dp0) - j = bc.dp0 - 1; + if (--j <= dp1 && j >= dp0) + j = dp0 - 1; if (s0[j] != '0') break; --i; @@ -1603,11 +1611,11 @@ _Py_dg_strtod(const char *s00, char **se) y = 0; for(i = 0; i < nd0; ++i) y = 10*y + s0[i] - '0'; - for(j = bc.dp1; i < nd; ++i) + for(j = dp1; i < nd; ++i) y = 10*y + s0[j++] - '0'; } } - bd0 = s2b(s0, nd0, nd, y, bc.dplen); + bd0 = s2b(s0, nd0, nd, y); if (bd0 == NULL) goto failed_malloc; @@ -1730,6 +1738,30 @@ _Py_dg_strtod(const char *s00, char **se) if (bc.nd > nd && i <= 0) { if (bc.dsign) break; /* Must use bigcomp(). */ + + /* Here rv overestimates the truncated decimal value by at most + 0.5 ulp(rv). Hence rv either overestimates the true decimal + value by <= 0.5 ulp(rv), or underestimates it by some small + amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of + the true decimal value, so it's possible to exit. + + Exception: if scaled rv is a normal exact power of 2, but not + DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the + next double, so the correctly rounded result is either rv - 0.5 + ulp(rv) or rv; in this case, use bigcomp to distinguish. */ + + if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) { + /* rv can't be 0, since it's an overestimate for some + nonzero value. So rv is a normal power of 2. */ + j = (int)(word0(&rv) & Exp_mask) >> Exp_shift; + /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if + rv / 2^bc.scale >= 2^-1021. */ + if (j - bc.scale >= 2) { + dval(&rv) -= 0.5 * sulp(&rv, &bc); + break; + } + } + { bc.nd = nd; i = -1; /* Discarded digits make delta smaller. */ -- cgit v0.12