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. ........
This commit is contained in:
parent
ae5465a578
commit
853c3bbc4c
|
@ -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<sign>[-+])? # an optional sign, followed by
|
||||||
|
(?=\d|\.\d) # a number with at least one digit
|
||||||
|
(?P<int>\d*) # having a (possibly empty) integer part
|
||||||
|
(?:\.(?P<frac>\d*))? # followed by an optional fractional part
|
||||||
|
(?:E(?P<exp>[-+]?\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('<d', struct.pack('<Q', bits + (bits&1)))[0]
|
||||||
|
#self.assertEqual(fs, float(s))
|
||||||
|
|
||||||
|
def test_boundaries(self):
|
||||||
|
# boundaries expressed as triples (n, e, u), where
|
||||||
|
# n*10**e is an approximation to the boundary value and
|
||||||
|
# u*10**e is 1ulp
|
||||||
|
boundaries = [
|
||||||
|
(10000000000000000000, -19, 1110), # a power of 2 boundary (1.0)
|
||||||
|
(17976931348623159077, 289, 1995), # overflow boundary (2.**1024)
|
||||||
|
(22250738585072013831, -327, 4941), # normal/subnormal (2.**-1022)
|
||||||
|
(0, -327, 4941), # zero
|
||||||
|
]
|
||||||
|
for n, e, u in boundaries:
|
||||||
|
for j in range(1000):
|
||||||
|
for i in range(TEST_SIZE):
|
||||||
|
digits = n + random.randrange(-3*u, 3*u)
|
||||||
|
exponent = e
|
||||||
|
s = '{}e{}'.format(digits, exponent)
|
||||||
|
self.check_strtod(s)
|
||||||
|
n *= 10
|
||||||
|
u *= 10
|
||||||
|
e -= 1
|
||||||
|
|
||||||
|
def test_underflow_boundary(self):
|
||||||
|
# test values close to 2**-1075, the underflow boundary; similar
|
||||||
|
# to boundary_tests, except that the random error doesn't scale
|
||||||
|
# with n
|
||||||
|
for exponent in range(-400, -320):
|
||||||
|
base = 10**-exponent // 2**1075
|
||||||
|
for j in range(TEST_SIZE):
|
||||||
|
digits = base + random.randrange(-1000, 1000)
|
||||||
|
s = '{}e{}'.format(digits, exponent)
|
||||||
|
self.check_strtod(s)
|
||||||
|
|
||||||
|
def test_bigcomp(self):
|
||||||
|
DIG10 = 10**50
|
||||||
|
for i in range(1000):
|
||||||
|
for j in range(TEST_SIZE):
|
||||||
|
digits = random.randrange(DIG10)
|
||||||
|
exponent = random.randrange(-400, 400)
|
||||||
|
s = '{}e{}'.format(digits, exponent)
|
||||||
|
self.check_strtod(s)
|
||||||
|
|
||||||
|
def test_parsing(self):
|
||||||
|
digits = tuple(map(str, range(10)))
|
||||||
|
signs = ('+', '-', '')
|
||||||
|
|
||||||
|
# put together random short valid strings
|
||||||
|
# \d*[.\d*]?e
|
||||||
|
for i in range(1000):
|
||||||
|
for j in range(TEST_SIZE):
|
||||||
|
s = random.choice(signs)
|
||||||
|
intpart_len = random.randrange(5)
|
||||||
|
s += ''.join(random.choice(digits) for _ in range(intpart_len))
|
||||||
|
if random.choice([True, False]):
|
||||||
|
s += '.'
|
||||||
|
fracpart_len = random.randrange(5)
|
||||||
|
s += ''.join(random.choice(digits)
|
||||||
|
for _ in range(fracpart_len))
|
||||||
|
else:
|
||||||
|
fracpart_len = 0
|
||||||
|
if random.choice([True, False]):
|
||||||
|
s += random.choice(['e', 'E'])
|
||||||
|
s += random.choice(signs)
|
||||||
|
exponent_len = random.randrange(1, 4)
|
||||||
|
s += ''.join(random.choice(digits)
|
||||||
|
for _ in range(exponent_len))
|
||||||
|
|
||||||
|
if intpart_len + fracpart_len:
|
||||||
|
self.check_strtod(s)
|
||||||
|
else:
|
||||||
|
try:
|
||||||
|
float(s)
|
||||||
|
except ValueError:
|
||||||
|
pass
|
||||||
|
else:
|
||||||
|
assert False, "expected ValueError"
|
||||||
|
|
||||||
|
def test_particular(self):
|
||||||
|
# inputs that produced crashes or incorrectly rounded results with
|
||||||
|
# previous versions of dtoa.c, for various reasons
|
||||||
|
test_strings = [
|
||||||
|
# issue 7632 bug 1, originally reported failing case
|
||||||
|
'2183167012312112312312.23538020374420446192e-370',
|
||||||
|
# 5 instances of issue 7632 bug 2
|
||||||
|
'12579816049008305546974391768996369464963024663104e-357',
|
||||||
|
'17489628565202117263145367596028389348922981857013e-357',
|
||||||
|
'18487398785991994634182916638542680759613590482273e-357',
|
||||||
|
'32002864200581033134358724675198044527469366773928e-358',
|
||||||
|
'94393431193180696942841837085033647913224148539854e-358',
|
||||||
|
# failing case for bug introduced by METD in r77451 (attempted
|
||||||
|
# fix for issue 7632, bug 2), and fixed in r77482.
|
||||||
|
'28639097178261763178489759107321392745108491825303e-311',
|
||||||
|
# two numbers demonstrating a flaw in the bigcomp 'dig == 0'
|
||||||
|
# correction block (issue 7632, bug 3)
|
||||||
|
'1.00000000000000001e44',
|
||||||
|
'1.0000000000000000100000000000000000000001e44',
|
||||||
|
# dtoa.c bug for numbers just smaller than a power of 2 (issue
|
||||||
|
# 7632, bug 4)
|
||||||
|
'99999999999999994487665465554760717039532578546e-47',
|
||||||
|
# failing case for off-by-one error introduced by METD in
|
||||||
|
# r77483 (dtoa.c cleanup), fixed in r77490
|
||||||
|
'965437176333654931799035513671997118345570045914469' #...
|
||||||
|
'6213413350821416312194420007991306908470147322020121018368e0',
|
||||||
|
# incorrect lsb detection for round-half-to-even when
|
||||||
|
# bc->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()
|
266
Python/dtoa.c
266
Python/dtoa.c
|
@ -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 dp0, dp1, dplen, dsign, e0, nd, nd0, scale;
|
int dsign, e0, nd, nd0, scale;
|
||||||
};
|
};
|
||||||
|
|
||||||
#define FFFFFFFF 0xffffffffUL
|
#define FFFFFFFF 0xffffffffUL
|
||||||
|
@ -437,7 +437,7 @@ multadd(Bigint *b, int m, int a) /* multiply by m and add a */
|
||||||
NULL on failure. */
|
NULL on failure. */
|
||||||
|
|
||||||
static Bigint *
|
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;
|
Bigint *b;
|
||||||
int i, k;
|
int i, k;
|
||||||
|
@ -451,18 +451,16 @@ s2b(const char *s, int nd0, int nd, ULong y9, int dplen)
|
||||||
b->x[0] = y9;
|
b->x[0] = y9;
|
||||||
b->wds = 1;
|
b->wds = 1;
|
||||||
|
|
||||||
i = 9;
|
if (nd <= 9)
|
||||||
if (9 < nd0) {
|
return b;
|
||||||
s += 9;
|
|
||||||
do {
|
s += 9;
|
||||||
b = multadd(b, 10, *s++ - '0');
|
for (i = 9; i < nd0; i++) {
|
||||||
if (b == NULL)
|
b = multadd(b, 10, *s++ - '0');
|
||||||
return NULL;
|
if (b == NULL)
|
||||||
} while(++i < nd0);
|
return NULL;
|
||||||
s += dplen;
|
|
||||||
}
|
}
|
||||||
else
|
s++;
|
||||||
s += dplen + 9;
|
|
||||||
for(; i < nd; i++) {
|
for(; i < nd; i++) {
|
||||||
b = multadd(b, 10, *s++ - '0');
|
b = multadd(b, 10, *s++ - '0');
|
||||||
if (b == NULL)
|
if (b == NULL)
|
||||||
|
@ -1130,76 +1128,120 @@ quorem(Bigint *b, Bigint *S)
|
||||||
return q;
|
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
|
Assuming that x is finite and nonnegative (positive zero is fine
|
||||||
representable as a double, sulp(x) is equivalent to 2^bc.scale * ulp(x /
|
here) and x / 2^bc.scale is exactly representable as a double,
|
||||||
2^bc.scale). */
|
sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */
|
||||||
|
|
||||||
static double
|
static double
|
||||||
sulp(U *x, BCinfo *bc)
|
sulp(U *x, BCinfo *bc)
|
||||||
{
|
{
|
||||||
U u;
|
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 */
|
/* rv/2^bc->scale is subnormal */
|
||||||
word0(&u) = (P+2)*Exp_msk1;
|
word0(&u) = (P+2)*Exp_msk1;
|
||||||
word1(&u) = 0;
|
word1(&u) = 0;
|
||||||
return u.d;
|
return u.d;
|
||||||
}
|
}
|
||||||
else
|
else {
|
||||||
|
assert(word0(x) || word1(x)); /* x != 0.0 */
|
||||||
return ulp(x);
|
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
|
static int
|
||||||
bigcomp(U *rv, const char *s0, BCinfo *bc)
|
bigcomp(U *rv, const char *s0, BCinfo *bc)
|
||||||
{
|
{
|
||||||
Bigint *b, *d;
|
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;
|
nd = bc->nd;
|
||||||
nd0 = bc->nd0;
|
nd0 = bc->nd0;
|
||||||
p5 = nd + bc->e0;
|
p5 = nd + bc->e0;
|
||||||
speccase = 0;
|
if (rv->d == 0.) {
|
||||||
if (rv->d == 0.) { /* special case: value near underflow-to-zero */
|
/* special case because d2b doesn't handle 0.0 */
|
||||||
/* threshold was rounded to zero */
|
b = i2b(0);
|
||||||
b = i2b(1);
|
|
||||||
if (b == NULL)
|
if (b == NULL)
|
||||||
return -1;
|
return -1;
|
||||||
p2 = Emin - P + 1;
|
p2 = Emin - P + 1; /* = -1074 for IEEE 754 binary64 */
|
||||||
bbits = 1;
|
bbits = 0;
|
||||||
word0(rv) = (P+2) << Exp_shift;
|
|
||||||
i = 0;
|
|
||||||
{
|
|
||||||
speccase = 1;
|
|
||||||
--p2;
|
|
||||||
dsign = 0;
|
|
||||||
goto have_i;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
else
|
else {
|
||||||
{
|
|
||||||
b = d2b(rv, &p2, &bbits);
|
b = d2b(rv, &p2, &bbits);
|
||||||
if (b == NULL)
|
if (b == NULL)
|
||||||
return -1;
|
return -1;
|
||||||
|
p2 -= bc->scale;
|
||||||
}
|
}
|
||||||
p2 -= bc->scale;
|
/* now rv/2^(bc->scale) = b * 2**p2, and b has bbits significant bits */
|
||||||
/* floor(log2(rv)) == bbits - 1 + p2 */
|
|
||||||
/* Check for denormal case. */
|
/* 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;
|
i = P - bbits;
|
||||||
if (i > (j = P - Emin - 1 + p2)) {
|
if (i > p2 - (Emin - P + 1))
|
||||||
i = j;
|
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 = lshift(b, ++i);
|
||||||
if (b == NULL)
|
if (b == NULL)
|
||||||
return -1;
|
return -1;
|
||||||
b->x[0] |= 1;
|
/* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway
|
||||||
}
|
case, this is used for round to even. */
|
||||||
have_i:
|
odd = b->x[0] & 2;
|
||||||
|
b->x[0] |= 1;
|
||||||
|
|
||||||
p2 -= p5 + i;
|
p2 -= p5 + i;
|
||||||
d = i2b(1);
|
d = i2b(1);
|
||||||
if (d == NULL) {
|
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
|
/* if b >= d, round down */
|
||||||
on either side of the input string. If b >= d, round down. */
|
|
||||||
if (cmp(b, d) >= 0) {
|
if (cmp(b, d) >= 0) {
|
||||||
dd = -1;
|
dd = -1;
|
||||||
goto ret;
|
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 */
|
/* Compare b/d with s0 */
|
||||||
|
for(i = 0; i < nd0; i++) {
|
||||||
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;
|
|
||||||
}
|
|
||||||
b = multadd(b, 10, 0);
|
b = multadd(b, 10, 0);
|
||||||
if (b == NULL) {
|
if (b == NULL) {
|
||||||
Bfree(d);
|
Bfree(d);
|
||||||
return -1;
|
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;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
for(j = bc->dp1; i++ < nd;) {
|
s0++;
|
||||||
if ((dd = s0[j++] - '0' - dig))
|
for(; i < nd; i++) {
|
||||||
goto ret;
|
|
||||||
if (!b->x[0] && b->wds == 1) {
|
|
||||||
if (i < nd)
|
|
||||||
dd = 1;
|
|
||||||
goto ret;
|
|
||||||
}
|
|
||||||
b = multadd(b, 10, 0);
|
b = multadd(b, 10, 0);
|
||||||
if (b == NULL) {
|
if (b == NULL) {
|
||||||
Bfree(d);
|
Bfree(d);
|
||||||
return -1;
|
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)
|
if (b->x[0] || b->wds > 1)
|
||||||
dd = -1;
|
dd = -1;
|
||||||
ret:
|
ret:
|
||||||
Bfree(b);
|
Bfree(b);
|
||||||
Bfree(d);
|
Bfree(d);
|
||||||
if (speccase) {
|
if (dd > 0 || (dd == 0 && odd))
|
||||||
if (dd <= 0)
|
dval(rv) += sulp(rv, bc);
|
||||||
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;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return 0;
|
return 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
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, dp0, dp1, dplen, e, e1, error;
|
||||||
int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
|
int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
|
||||||
const char *s, *s0, *s1;
|
const char *s, *s0, *s1;
|
||||||
double aadj, aadj1;
|
double aadj, aadj1;
|
||||||
|
@ -1341,7 +1349,7 @@ _Py_dg_strtod(const char *s00, char **se)
|
||||||
BCinfo bc;
|
BCinfo bc;
|
||||||
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
|
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
|
||||||
|
|
||||||
sign = nz0 = nz = bc.dplen = 0;
|
sign = nz0 = nz = dplen = 0;
|
||||||
dval(&rv) = 0.;
|
dval(&rv) = 0.;
|
||||||
for(s = s00;;s++) switch(*s) {
|
for(s = s00;;s++) switch(*s) {
|
||||||
case '-':
|
case '-':
|
||||||
|
@ -1380,11 +1388,11 @@ _Py_dg_strtod(const char *s00, char **se)
|
||||||
else if (nd < 16)
|
else if (nd < 16)
|
||||||
z = 10*z + c - '0';
|
z = 10*z + c - '0';
|
||||||
nd0 = nd;
|
nd0 = nd;
|
||||||
bc.dp0 = bc.dp1 = s - s0;
|
dp0 = dp1 = s - s0;
|
||||||
if (c == '.') {
|
if (c == '.') {
|
||||||
c = *++s;
|
c = *++s;
|
||||||
bc.dp1 = s - s0;
|
dp1 = s - s0;
|
||||||
bc.dplen = bc.dp1 - bc.dp0;
|
dplen = 1;
|
||||||
if (!nd) {
|
if (!nd) {
|
||||||
for(; c == '0'; c = *++s)
|
for(; c == '0'; c = *++s)
|
||||||
nz++;
|
nz++;
|
||||||
|
@ -1587,10 +1595,10 @@ _Py_dg_strtod(const char *s00, char **se)
|
||||||
/* in IEEE arithmetic. */
|
/* in IEEE arithmetic. */
|
||||||
i = j = 18;
|
i = j = 18;
|
||||||
if (i > nd0)
|
if (i > nd0)
|
||||||
j += bc.dplen;
|
j += dplen;
|
||||||
for(;;) {
|
for(;;) {
|
||||||
if (--j <= bc.dp1 && j >= bc.dp0)
|
if (--j <= dp1 && j >= dp0)
|
||||||
j = bc.dp0 - 1;
|
j = dp0 - 1;
|
||||||
if (s0[j] != '0')
|
if (s0[j] != '0')
|
||||||
break;
|
break;
|
||||||
--i;
|
--i;
|
||||||
|
@ -1603,11 +1611,11 @@ _Py_dg_strtod(const char *s00, char **se)
|
||||||
y = 0;
|
y = 0;
|
||||||
for(i = 0; i < nd0; ++i)
|
for(i = 0; i < nd0; ++i)
|
||||||
y = 10*y + s0[i] - '0';
|
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';
|
y = 10*y + s0[j++] - '0';
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
bd0 = s2b(s0, nd0, nd, y, bc.dplen);
|
bd0 = s2b(s0, nd0, nd, y);
|
||||||
if (bd0 == NULL)
|
if (bd0 == NULL)
|
||||||
goto failed_malloc;
|
goto failed_malloc;
|
||||||
|
|
||||||
|
@ -1730,6 +1738,30 @@ _Py_dg_strtod(const char *s00, char **se)
|
||||||
if (bc.nd > nd && i <= 0) {
|
if (bc.nd > nd && i <= 0) {
|
||||||
if (bc.dsign)
|
if (bc.dsign)
|
||||||
break; /* Must use bigcomp(). */
|
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;
|
bc.nd = nd;
|
||||||
i = -1; /* Discarded digits make delta smaller. */
|
i = -1; /* Discarded digits make delta smaller. */
|
||||||
|
|
Loading…
Reference in New Issue