Author: mark.dickinson Date: Sun Jan 24 11:16:29 2010 New Revision: 77722 Log: Merged revisions 77691,77698,77713-77714 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r77691 | mark.dickinson | 2010-01-22 16:18:09 +0000 (Fri, 22 Jan 2010) | 1 line Correct typo in comment. ........ r77698 | mark.dickinson | 2010-01-22 17:04:07 +0000 (Fri, 22 Jan 2010) | 3 lines Issue #7743: Fix a potential incorrect rounding bug in dtoa.c (2nd bug in issue 7743). ........ r77713 | mark.dickinson | 2010-01-23 20:48:56 +0000 (Sat, 23 Jan 2010) | 3 lines Issue #7743: Add checks for zero inputs to the lshift and mult functions; this fixes the first bug described in issue #7743. ........ r77714 | mark.dickinson | 2010-01-23 21:25:53 +0000 (Sat, 23 Jan 2010) | 1 line dtoa.c fix from upstream that fixes incorrectly rounded results for certain subnormals that are also halfway cases. ........ Modified: python/branches/py3k/ (props changed) python/branches/py3k/Python/dtoa.c python/branches/py3k/Python/pystrtod.c Modified: python/branches/py3k/Python/dtoa.c ============================================================================== --- python/branches/py3k/Python/dtoa.c (original) +++ python/branches/py3k/Python/dtoa.c Sun Jan 24 11:16:29 2010 @@ -235,6 +235,7 @@ #define Bias 1023 #define Emax 1023 #define Emin (-1022) +#define Etiny (-1074) /* smallest denormal is 2**Etiny */ #define Exp_1 0x3ff00000 #define Exp_11 0x3ff00000 #define Ebits 11 @@ -244,7 +245,6 @@ #define Bletch 0x10 #define Bndry_mask 0xfffff #define Bndry_mask1 0xfffff -#define LSB 1 #define Sign_bit 0x80000000 #define Log2P 1 #define Tiny0 0 @@ -622,6 +622,15 @@ ULong z2; #endif + if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) { + c = Balloc(0); + if (c == NULL) + return NULL; + c->wds = 1; + c->x[0] = 0; + return c; + } + if (a->wds < b->wds) { c = a; a = b; @@ -820,6 +829,9 @@ Bigint *b1; ULong *x, *x1, *xe, z; + if (!k || (!b->x[0] && b->wds == 1)) + return b; + n = k >> 5; k1 = b->k; n1 = n + b->wds + 1; @@ -1019,6 +1031,76 @@ return dval(&d); } +/* Convert a scaled double to a Bigint plus an exponent. Similar to d2b, + except that it accepts the scale parameter used in _Py_dg_strtod (which + should be either 0 or 2*P), and the normalization for the return value is + different (see below). On input, d should be finite and nonnegative, and d + / 2**scale should be exactly representable as an IEEE 754 double. + + Returns a Bigint b and an integer e such that + + dval(d) / 2**scale = b * 2**e. + + Unlike d2b, b is not necessarily odd: b and e are normalized so + that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P + and e == Etiny. This applies equally to an input of 0.0: in that + case the return values are b = 0 and e = Etiny. + + The above normalization ensures that for all possible inputs d, + 2**e gives ulp(d/2**scale). + + Returns NULL on failure. +*/ + +static Bigint * +sd2b(U *d, int scale, int *e) +{ + Bigint *b; + + b = Balloc(1); + if (b == NULL) + return NULL; + + /* First construct b and e assuming that scale == 0. */ + b->wds = 2; + b->x[0] = word1(d); + b->x[1] = word0(d) & Frac_mask; + *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift); + if (*e < Etiny) + *e = Etiny; + else + b->x[1] |= Exp_msk1; + + /* Now adjust for scale, provided that b != 0. */ + if (scale && (b->x[0] || b->x[1])) { + *e -= scale; + if (*e < Etiny) { + scale = Etiny - *e; + *e = Etiny; + /* We can't shift more than P-1 bits without shifting out a 1. */ + assert(0 < scale && scale <= P - 1); + if (scale >= 32) { + /* The bits shifted out should all be zero. */ + assert(b->x[0] == 0); + b->x[0] = b->x[1]; + b->x[1] = 0; + scale -= 32; + } + if (scale) { + /* The bits shifted out should all be zero. */ + assert(b->x[0] << (32 - scale) == 0); + b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale)); + b->x[1] >>= scale; + } + } + } + /* Ensure b is normalized. */ + if (!b->x[1]) + b->wds = 1; + + return b; +} + /* Convert a double to a Bigint plus an exponent. Return NULL on failure. Given a finite nonzero double d, return an odd Bigint b and exponent *e @@ -1028,7 +1110,6 @@ If d is zero, then b == 0, *e == -1010, *bbits = 0. */ - static Bigint * d2b(U *d, int *e, int *bits) { @@ -1299,45 +1380,29 @@ bigcomp(U *rv, const char *s0, BCinfo *bc) { Bigint *b, *d; - int b2, bbits, d2, dd, i, nd, nd0, odd, p2, p5; + int b2, d2, dd, i, nd, nd0, odd, p2, p5; dd = 0; /* silence compiler warning about possibly unused variable */ nd = bc->nd; nd0 = bc->nd0; p5 = nd + bc->e0; - 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; /* = -1074 for IEEE 754 binary64 */ - bbits = 0; - } - else { - b = d2b(rv, &p2, &bbits); - if (b == NULL) - return -1; - p2 -= bc->scale; - } - /* 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 > 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); + b = sd2b(rv, bc->scale, &p2); 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; + odd = b->x[0] & 1; + + /* left shift b by 1 bit and or a 1 into the least significant bit; + this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */ + b = lshift(b, 1); + if (b == NULL) + return -1; b->x[0] |= 1; + p2--; - p2 -= p5 + i; + p2 -= p5; d = i2b(1); if (d == NULL) { Bfree(b); @@ -1425,8 +1490,8 @@ double _Py_dg_strtod(const char *s00, char **se) { - int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e, e1, error; - int esign, i, j, k, lz, nd, nd0, sign; + int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error; + int esign, i, j, k, lz, nd, nd0, odd, sign; const char *s, *s0, *s1; double aadj, aadj1; U aadj2, adj, rv, rv0; @@ -1786,13 +1851,17 @@ goto failed_malloc; } Bcopy(bd, bd0); - bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ + bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */ if (bb == NULL) { Bfree(bd); Bfree(bd0); goto failed_malloc; } - /* tdv = bd * 10^e; srv = bb * 2^(bbe - scale) */ + /* Record whether lsb of bb is odd, in case we need this + for the round-to-even step later. */ + odd = bb->x[0] & 1; + + /* tdv = bd * 10**e; srv = bb * 2**bbe */ bs = i2b(1); if (bs == NULL) { Bfree(bb); @@ -1813,44 +1882,28 @@ bb2 += bbe; else bd2 -= bbe; + bs2 = bb2; + bb2++; + bd2++; - /* At this stage e = bd2 - bb2 + bbe = bd5 - bb5, so: - - tdv = bd * 2^(bbe + bd2 - bb2) * 5^(bd5 - bb5) - srv = bb * 2^(bbe - scale). + /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1, + and bs == 1, so: - Compute j such that + tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5) + srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2) + 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2) - 0.5 ulp(srv) = 2^(bbe - scale - j) - */ + It follows that: - bs2 = bb2; - j = bbe - bc.scale; - i = j + bbbits - 1; /* logb(rv) */ - if (i < Emin) /* denormal */ - j += P - Emin; - else - j = P + 1 - bbbits; + M * tdv = bd * 2**bd2 * 5**bd5 + M * srv = bb * 2**bb2 * 5**bb5 + M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5 - /* 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. + for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but + this fact is not needed below.) */ - bb2 += j; - bd2 += j; - 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). */ + /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */ i = bb2 < bd2 ? bb2 : bd2; if (i > bs2) i = bs2; @@ -2028,12 +2081,12 @@ word1(&rv) = 0xffffffff; break; } - if (!(word1(&rv) & LSB)) + if (!odd) break; if (dsign) - dval(&rv) += ulp(&rv); + dval(&rv) += sulp(&rv, &bc); else { - dval(&rv) -= ulp(&rv); + dval(&rv) -= sulp(&rv, &bc); if (!dval(&rv)) { if (bc.nd >nd) break; Modified: python/branches/py3k/Python/pystrtod.c ============================================================================== --- python/branches/py3k/Python/pystrtod.c (original) +++ python/branches/py3k/Python/pystrtod.c Sun Jan 24 11:16:29 2010 @@ -104,7 +104,7 @@ _Py_SET_53BIT_PRECISION_END; if (*endptr == nptr) - /* string might represent and inf or nan */ + /* string might represent an inf or nan */ result = _Py_parse_inf_or_nan(nptr, endptr); return result;