diff --git a/Lib/test/ieee754.txt b/Lib/test/ieee754.txt new file mode 100644 index 0000000000..3e986cdb10 --- /dev/null +++ b/Lib/test/ieee754.txt @@ -0,0 +1,183 @@ +====================================== +Python IEEE 754 floating point support +====================================== + +>>> from sys import float_info as FI +>>> from math import * +>>> PI = pi +>>> E = e + +You must never compare two floats with == because you are not going to get +what you expect. We treat two floats as equal if the difference between them +is small than epsilon. +>>> EPS = 1E-15 +>>> def equal(x, y): +... """Almost equal helper for floats""" +... return abs(x - y) < EPS + + +NaNs and INFs +============= + +In Python 2.6 and newer NaNs (not a number) and infinity can be constructed +from the strings 'inf' and 'nan'. + +>>> INF = float('inf') +>>> NINF = float('-inf') +>>> NAN = float('nan') + +>>> INF +inf +>>> NINF +-inf +>>> NAN +nan + +The math module's ``isnan`` and ``isinf`` functions can be used to detect INF +and NAN: +>>> isinf(INF), isinf(NINF), isnan(NAN) +(True, True, True) +>>> INF == -NINF +True + +Infinity +-------- + +Ambiguous operations like ``0 * inf`` or ``inf - inf`` result in NaN. +>>> INF * 0 +nan +>>> INF - INF +nan +>>> INF / INF +nan + +However unambiguous operations with inf return inf: +>>> INF * INF +inf +>>> 1.5 * INF +inf +>>> 0.5 * INF +inf +>>> INF / 1000 +inf + +Not a Number +------------ + +NaNs are never equal to another number, even itself +>>> NAN == NAN +False +>>> NAN < 0 +False +>>> NAN >= 0 +False + +All operations involving a NaN return a NaN except for nan**0 and 1**nan. +>>> 1 + NAN +nan +>>> 1 * NAN +nan +>>> 0 * NAN +nan +>>> 1 ** NAN +1.0 +>>> NAN ** 0 +1.0 +>>> 0 ** NAN +nan +>>> (1.0 + FI.epsilon) * NAN +nan + +Misc Functions +============== + +The power of 1 raised to x is always 1.0, even for special values like 0, +infinity and NaN. + +>>> pow(1, 0) +1.0 +>>> pow(1, INF) +1.0 +>>> pow(1, -INF) +1.0 +>>> pow(1, NAN) +1.0 + +The power of 0 raised to x is defined as 0, if x is positive. Negative +finite values are a domain error or zero division error and NaN result in a +silent NaN. + +>>> pow(0, 0) +1.0 +>>> pow(0, INF) +0.0 +>>> pow(0, -INF) +inf +>>> 0 ** -1 +Traceback (most recent call last): +... +ZeroDivisionError: 0.0 cannot be raised to a negative power +>>> pow(0, NAN) +nan + + +Trigonometric Functions +======================= + +>>> sin(INF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> sin(NINF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> sin(NAN) +nan +>>> cos(INF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> cos(NINF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> cos(NAN) +nan +>>> tan(INF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> tan(NINF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> tan(NAN) +nan + +Neither pi nor tan are exact, but you can assume that tan(pi/2) is a large value +and tan(pi) is a very small value: +>>> tan(PI/2) > 1E10 +True +>>> -tan(-PI/2) > 1E10 +True +>>> tan(PI) < 1E-15 +True + +>>> asin(NAN), acos(NAN), atan(NAN) +(nan, nan, nan) +>>> asin(INF), asin(NINF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> acos(INF), acos(NINF) +Traceback (most recent call last): +... +ValueError: math domain error +>>> equal(atan(INF), PI/2), equal(atan(NINF), -PI/2) +(True, True) + + +Hyberbolic Functions +==================== + diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 2f64180652..fa79456ed4 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -4,6 +4,7 @@ from test.support import verbose, requires_IEEE_754 from test import support import unittest +import fractions import itertools import decimal import math @@ -186,6 +187,9 @@ def result_check(expected, got, ulp_tol=5, abs_tol=0.0): # Check exactly equal (applies also to strings representing exceptions) if got == expected: + if not got and not expected: + if math.copysign(1, got) != math.copysign(1, expected): + return f"expected {expected}, got {got} (zero has wrong sign)" return None failure = "not equal" @@ -234,6 +238,10 @@ def __init__(self, value): def __index__(self): return self.value +class BadDescr: + def __get__(self, obj, objtype=None): + raise ValueError + class MathTests(unittest.TestCase): def ftest(self, name, got, expected, ulp_tol=5, abs_tol=0.0): @@ -323,6 +331,7 @@ def testAtan2(self): self.ftest('atan2(0, 1)', math.atan2(0, 1), 0) self.ftest('atan2(1, 1)', math.atan2(1, 1), math.pi/4) self.ftest('atan2(1, 0)', math.atan2(1, 0), math.pi/2) + self.ftest('atan2(1, -1)', math.atan2(1, -1), 3*math.pi/4) # math.atan2(0, x) self.ftest('atan2(0., -inf)', math.atan2(0., NINF), math.pi) @@ -416,16 +425,22 @@ def __ceil__(self): return 42 class TestNoCeil: pass + class TestBadCeil: + __ceil__ = BadDescr() self.assertEqual(math.ceil(TestCeil()), 42) self.assertEqual(math.ceil(FloatCeil()), 42) self.assertEqual(math.ceil(FloatLike(42.5)), 43) self.assertRaises(TypeError, math.ceil, TestNoCeil()) + self.assertRaises(ValueError, math.ceil, TestBadCeil()) t = TestNoCeil() t.__ceil__ = lambda *args: args self.assertRaises(TypeError, math.ceil, t) self.assertRaises(TypeError, math.ceil, t, 0) + self.assertEqual(math.ceil(FloatLike(+1.0)), +1.0) + self.assertEqual(math.ceil(FloatLike(-1.0)), -1.0) + @requires_IEEE_754 def testCopysign(self): self.assertEqual(math.copysign(1, 42), 1.0) @@ -566,16 +581,22 @@ def __floor__(self): return 42 class TestNoFloor: pass + class TestBadFloor: + __floor__ = BadDescr() self.assertEqual(math.floor(TestFloor()), 42) self.assertEqual(math.floor(FloatFloor()), 42) self.assertEqual(math.floor(FloatLike(41.9)), 41) self.assertRaises(TypeError, math.floor, TestNoFloor()) + self.assertRaises(ValueError, math.floor, TestBadFloor()) t = TestNoFloor() t.__floor__ = lambda *args: args self.assertRaises(TypeError, math.floor, t) self.assertRaises(TypeError, math.floor, t, 0) + self.assertEqual(math.floor(FloatLike(+1.0)), +1.0) + self.assertEqual(math.floor(FloatLike(-1.0)), -1.0) + def testFmod(self): self.assertRaises(TypeError, math.fmod) self.ftest('fmod(10, 1)', math.fmod(10, 1), 0.0) @@ -597,6 +618,7 @@ def testFmod(self): self.assertEqual(math.fmod(-3.0, NINF), -3.0) self.assertEqual(math.fmod(0.0, 3.0), 0.0) self.assertEqual(math.fmod(0.0, NINF), 0.0) + self.assertRaises(ValueError, math.fmod, INF, INF) def testFrexp(self): self.assertRaises(TypeError, math.frexp) @@ -638,7 +660,7 @@ def testFsum(self): def msum(iterable): """Full precision summation. Compute sum(iterable) without any intermediate accumulation of error. Based on the 'lsum' function - at http://code.activestate.com/recipes/393090/ + at https://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/ """ tmant, texp = 0, 0 @@ -666,6 +688,7 @@ def msum(iterable): ([], 0.0), ([0.0], 0.0), ([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100), + ([1e100, 1.0, -1e100, 1e-100, 1e50, -1, -1e50], 1e-100), ([2.0**53, -0.5, -2.0**-54], 2.0**53-1.0), ([2.0**53, 1.0, 2.0**-100], 2.0**53+2.0), ([2.0**53+10.0, 1.0, 2.0**-100], 2.0**53+12.0), @@ -713,6 +736,22 @@ def msum(iterable): s = msum(vals) self.assertEqual(msum(vals), math.fsum(vals)) + self.assertEqual(math.fsum([1.0, math.inf]), math.inf) + self.assertTrue(math.isnan(math.fsum([math.nan, 1.0]))) + self.assertEqual(math.fsum([1e100, FloatLike(1.0), -1e100, 1e-100, + 1e50, FloatLike(-1.0), -1e50]), 1e-100) + self.assertRaises(OverflowError, math.fsum, [1e+308, 1e+308]) + self.assertRaises(ValueError, math.fsum, [math.inf, -math.inf]) + self.assertRaises(TypeError, math.fsum, ['spam']) + self.assertRaises(TypeError, math.fsum, 1) + self.assertRaises(OverflowError, math.fsum, [10**1000]) + + def bad_iter(): + yield 1.0 + raise ZeroDivisionError + + self.assertRaises(ZeroDivisionError, math.fsum, bad_iter()) + def testGcd(self): gcd = math.gcd self.assertEqual(gcd(0, 0), 0) @@ -773,9 +812,13 @@ def testHypot(self): # Test allowable types (those with __float__) self.assertEqual(hypot(12.0, 5.0), 13.0) self.assertEqual(hypot(12, 5), 13) + self.assertEqual(hypot(0.75, -1), 1.25) + self.assertEqual(hypot(-1, 0.75), 1.25) + self.assertEqual(hypot(0.75, FloatLike(-1.)), 1.25) + self.assertEqual(hypot(FloatLike(-1.), 0.75), 1.25) self.assertEqual(hypot(Decimal(12), Decimal(5)), 13) self.assertEqual(hypot(Fraction(12, 32), Fraction(5, 32)), Fraction(13, 32)) - self.assertEqual(hypot(bool(1), bool(0), bool(1), bool(1)), math.sqrt(3)) + self.assertEqual(hypot(True, False, True, True, True), 2.0) # Test corner cases self.assertEqual(hypot(0.0, 0.0), 0.0) # Max input is zero @@ -830,6 +873,8 @@ def testHypot(self): scale = FLOAT_MIN / 2.0 ** exp self.assertEqual(math.hypot(4*scale, 3*scale), 5*scale) + self.assertRaises(TypeError, math.hypot, *([1.0]*18), 'spam') + @requires_IEEE_754 @unittest.skipIf(HAVE_DOUBLE_ROUNDING, "hypot() loses accuracy on machines with double rounding") @@ -922,12 +967,16 @@ def testDist(self): # Test allowable types (those with __float__) self.assertEqual(dist((14.0, 1.0), (2.0, -4.0)), 13.0) self.assertEqual(dist((14, 1), (2, -4)), 13) + self.assertEqual(dist((FloatLike(14.), 1), (2, -4)), 13) + self.assertEqual(dist((11, 1), (FloatLike(-1.), -4)), 13) + self.assertEqual(dist((14, FloatLike(-1.)), (2, -6)), 13) + self.assertEqual(dist((14, -1), (2, -6)), 13) self.assertEqual(dist((D(14), D(1)), (D(2), D(-4))), D(13)) self.assertEqual(dist((F(14, 32), F(1, 32)), (F(2, 32), F(-4, 32))), F(13, 32)) - self.assertEqual(dist((True, True, False, True, False), - (True, False, True, True, False)), - sqrt(2.0)) + self.assertEqual(dist((True, True, False, False, True, True), + (True, False, True, False, False, False)), + 2.0) # Test corner cases self.assertEqual(dist((13.25, 12.5, -3.25), @@ -965,6 +1014,8 @@ class T(tuple): dist((1, 2, 3, 4), (5, 6, 7)) with self.assertRaises(ValueError): # Check dimension agree dist((1, 2, 3), (4, 5, 6, 7)) + with self.assertRaises(TypeError): + dist((1,)*17 + ("spam",), (1,)*18) with self.assertRaises(TypeError): # Rejects invalid types dist("abc", "xyz") int_too_big_for_float = 10 ** (sys.float_info.max_10_exp + 5) @@ -972,6 +1023,16 @@ class T(tuple): dist((1, int_too_big_for_float), (2, 3)) with self.assertRaises((ValueError, OverflowError)): dist((2, 3), (1, int_too_big_for_float)) + with self.assertRaises(TypeError): + dist((1,), 2) + with self.assertRaises(TypeError): + dist([1], 2) + + class BadFloat: + __float__ = BadDescr() + + with self.assertRaises(ValueError): + dist([1], [BadFloat()]) # Verify that the one dimensional case is equivalent to abs() for i in range(20): @@ -1110,6 +1171,7 @@ def test_lcm(self): def testLdexp(self): self.assertRaises(TypeError, math.ldexp) + self.assertRaises(TypeError, math.ldexp, 2.0, 1.1) self.ftest('ldexp(0,1)', math.ldexp(0,1), 0) self.ftest('ldexp(1,1)', math.ldexp(1,1), 2) self.ftest('ldexp(1,-1)', math.ldexp(1,-1), 0.5) @@ -1142,6 +1204,7 @@ def testLdexp(self): def testLog(self): self.assertRaises(TypeError, math.log) + self.assertRaises(TypeError, math.log, 1, 2, 3) self.ftest('log(1/e)', math.log(1/math.e), -1) self.ftest('log(1)', math.log(1), 0) self.ftest('log(e)', math.log(math.e), 1) @@ -1152,6 +1215,7 @@ def testLog(self): 2302.5850929940457) self.assertRaises(ValueError, math.log, -1.5) self.assertRaises(ValueError, math.log, -10**1000) + self.assertRaises(ValueError, math.log, 10, -10) self.assertRaises(ValueError, math.log, NINF) self.assertEqual(math.log(INF), INF) self.assertTrue(math.isnan(math.log(NAN))) @@ -1202,6 +1266,277 @@ def testLog10(self): self.assertEqual(math.log(INF), INF) self.assertTrue(math.isnan(math.log10(NAN))) + def testSumProd(self): + sumprod = math.sumprod + Decimal = decimal.Decimal + Fraction = fractions.Fraction + + # Core functionality + self.assertEqual(sumprod(iter([10, 20, 30]), (1, 2, 3)), 140) + self.assertEqual(sumprod([1.5, 2.5], [3.5, 4.5]), 16.5) + self.assertEqual(sumprod([], []), 0) + self.assertEqual(sumprod([-1], [1.]), -1) + self.assertEqual(sumprod([1.], [-1]), -1) + + # Type preservation and coercion + for v in [ + (10, 20, 30), + (1.5, -2.5), + (Fraction(3, 5), Fraction(4, 5)), + (Decimal(3.5), Decimal(4.5)), + (2.5, 10), # float/int + (2.5, Fraction(3, 5)), # float/fraction + (25, Fraction(3, 5)), # int/fraction + (25, Decimal(4.5)), # int/decimal + ]: + for p, q in [(v, v), (v, v[::-1])]: + with self.subTest(p=p, q=q): + expected = sum(p_i * q_i for p_i, q_i in zip(p, q, strict=True)) + actual = sumprod(p, q) + self.assertEqual(expected, actual) + self.assertEqual(type(expected), type(actual)) + + # Bad arguments + self.assertRaises(TypeError, sumprod) # No args + self.assertRaises(TypeError, sumprod, []) # One arg + self.assertRaises(TypeError, sumprod, [], [], []) # Three args + self.assertRaises(TypeError, sumprod, None, [10]) # Non-iterable + self.assertRaises(TypeError, sumprod, [10], None) # Non-iterable + self.assertRaises(TypeError, sumprod, ['x'], [1.0]) + + # Uneven lengths + self.assertRaises(ValueError, sumprod, [10, 20], [30]) + self.assertRaises(ValueError, sumprod, [10], [20, 30]) + + # Overflows + self.assertEqual(sumprod([10**20], [1]), 10**20) + self.assertEqual(sumprod([1], [10**20]), 10**20) + self.assertEqual(sumprod([10**10], [10**10]), 10**20) + self.assertEqual(sumprod([10**7]*10**5, [10**7]*10**5), 10**19) + self.assertRaises(OverflowError, sumprod, [10**1000], [1.0]) + self.assertRaises(OverflowError, sumprod, [1.0], [10**1000]) + + # Error in iterator + def raise_after(n): + for i in range(n): + yield i + raise RuntimeError + with self.assertRaises(RuntimeError): + sumprod(range(10), raise_after(5)) + with self.assertRaises(RuntimeError): + sumprod(raise_after(5), range(10)) + + from test.test_iter import BasicIterClass + + self.assertEqual(sumprod(BasicIterClass(1), [1]), 0) + self.assertEqual(sumprod([1], BasicIterClass(1)), 0) + + # Error in multiplication + class BadMultiply: + def __mul__(self, other): + raise RuntimeError + def __rmul__(self, other): + raise RuntimeError + with self.assertRaises(RuntimeError): + sumprod([10, BadMultiply(), 30], [1, 2, 3]) + with self.assertRaises(RuntimeError): + sumprod([1, 2, 3], [10, BadMultiply(), 30]) + + # Error in addition + with self.assertRaises(TypeError): + sumprod(['abc', 3], [5, 10]) + with self.assertRaises(TypeError): + sumprod([5, 10], ['abc', 3]) + + # Special values should give the same as the pure python recipe + self.assertEqual(sumprod([10.1, math.inf], [20.2, 30.3]), math.inf) + self.assertEqual(sumprod([10.1, math.inf], [math.inf, 30.3]), math.inf) + self.assertEqual(sumprod([10.1, math.inf], [math.inf, math.inf]), math.inf) + self.assertEqual(sumprod([10.1, -math.inf], [20.2, 30.3]), -math.inf) + self.assertTrue(math.isnan(sumprod([10.1, math.inf], [-math.inf, math.inf]))) + self.assertTrue(math.isnan(sumprod([10.1, math.nan], [20.2, 30.3]))) + self.assertTrue(math.isnan(sumprod([10.1, math.inf], [math.nan, 30.3]))) + self.assertTrue(math.isnan(sumprod([10.1, math.inf], [20.3, math.nan]))) + + # Error cases that arose during development + args = ((-5, -5, 10), (1.5, 4611686018427387904, 2305843009213693952)) + self.assertEqual(sumprod(*args), 0.0) + + + @requires_IEEE_754 + @unittest.skipIf(HAVE_DOUBLE_ROUNDING, + "sumprod() accuracy not guaranteed on machines with double rounding") + @support.cpython_only # Other implementations may choose a different algorithm + def test_sumprod_accuracy(self): + sumprod = math.sumprod + self.assertEqual(sumprod([0.1] * 10, [1]*10), 1.0) + self.assertEqual(sumprod([0.1] * 20, [True, False] * 10), 1.0) + self.assertEqual(sumprod([True, False] * 10, [0.1] * 20), 1.0) + self.assertEqual(sumprod([1.0, 10E100, 1.0, -10E100], [1.0]*4), 2.0) + + @support.requires_resource('cpu') + def test_sumprod_stress(self): + sumprod = math.sumprod + product = itertools.product + Decimal = decimal.Decimal + Fraction = fractions.Fraction + + class Int(int): + def __add__(self, other): + return Int(int(self) + int(other)) + def __mul__(self, other): + return Int(int(self) * int(other)) + __radd__ = __add__ + __rmul__ = __mul__ + def __repr__(self): + return f'Int({int(self)})' + + class Flt(float): + def __add__(self, other): + return Int(int(self) + int(other)) + def __mul__(self, other): + return Int(int(self) * int(other)) + __radd__ = __add__ + __rmul__ = __mul__ + def __repr__(self): + return f'Flt({int(self)})' + + def baseline_sumprod(p, q): + """This defines the target behavior including exceptions and special values. + However, it is subject to rounding errors, so float inputs should be exactly + representable with only a few bits. + """ + total = 0 + for p_i, q_i in zip(p, q, strict=True): + total += p_i * q_i + return total + + def run(func, *args): + "Make comparing functions easier. Returns error status, type, and result." + try: + result = func(*args) + except (AssertionError, NameError): + raise + except Exception as e: + return type(e), None, 'None' + return None, type(result), repr(result) + + pools = [ + (-5, 10, -2**20, 2**31, 2**40, 2**61, 2**62, 2**80, 1.5, Int(7)), + (5.25, -3.5, 4.75, 11.25, 400.5, 0.046875, 0.25, -1.0, -0.078125), + (-19.0*2**500, 11*2**1000, -3*2**1500, 17*2*333, + 5.25, -3.25, -3.0*2**(-333), 3, 2**513), + (3.75, 2.5, -1.5, float('inf'), -float('inf'), float('NaN'), 14, + 9, 3+4j, Flt(13), 0.0), + (13.25, -4.25, Decimal('10.5'), Decimal('-2.25'), Fraction(13, 8), + Fraction(-11, 16), 4.75 + 0.125j, 97, -41, Int(3)), + (Decimal('6.125'), Decimal('12.375'), Decimal('-2.75'), Decimal(0), + Decimal('Inf'), -Decimal('Inf'), Decimal('NaN'), 12, 13.5), + (-2.0 ** -1000, 11*2**1000, 3, 7, -37*2**32, -2*2**-537, -2*2**-538, + 2*2**-513), + (-7 * 2.0 ** -510, 5 * 2.0 ** -520, 17, -19.0, -6.25), + (11.25, -3.75, -0.625, 23.375, True, False, 7, Int(5)), + ] + + for pool in pools: + for size in range(4): + for args1 in product(pool, repeat=size): + for args2 in product(pool, repeat=size): + args = (args1, args2) + self.assertEqual( + run(baseline_sumprod, *args), + run(sumprod, *args), + args, + ) + + @requires_IEEE_754 + @unittest.skipIf(HAVE_DOUBLE_ROUNDING, + "sumprod() accuracy not guaranteed on machines with double rounding") + @support.cpython_only # Other implementations may choose a different algorithm + @support.requires_resource('cpu') + def test_sumprod_extended_precision_accuracy(self): + import operator + from fractions import Fraction + from itertools import starmap + from collections import namedtuple + from math import log2, exp2, fabs + from random import choices, uniform, shuffle + from statistics import median + + DotExample = namedtuple('DotExample', ('x', 'y', 'target_sumprod', 'condition')) + + def DotExact(x, y): + vec1 = map(Fraction, x) + vec2 = map(Fraction, y) + return sum(starmap(operator.mul, zip(vec1, vec2, strict=True))) + + def Condition(x, y): + return 2.0 * DotExact(map(abs, x), map(abs, y)) / abs(DotExact(x, y)) + + def linspace(lo, hi, n): + width = (hi - lo) / (n - 1) + return [lo + width * i for i in range(n)] + + def GenDot(n, c): + """ Algorithm 6.1 (GenDot) works as follows. The condition number (5.7) of + the dot product xT y is proportional to the degree of cancellation. In + order to achieve a prescribed cancellation, we generate the first half of + the vectors x and y randomly within a large exponent range. This range is + chosen according to the anticipated condition number. The second half of x + and y is then constructed choosing xi randomly with decreasing exponent, + and calculating yi such that some cancellation occurs. Finally, we permute + the vectors x, y randomly and calculate the achieved condition number. + """ + + assert n >= 6 + n2 = n // 2 + x = [0.0] * n + y = [0.0] * n + b = log2(c) + + # First half with exponents from 0 to |_b/2_| and random ints in between + e = choices(range(int(b/2)), k=n2) + e[0] = int(b / 2) + 1 + e[-1] = 0.0 + + x[:n2] = [uniform(-1.0, 1.0) * exp2(p) for p in e] + y[:n2] = [uniform(-1.0, 1.0) * exp2(p) for p in e] + + # Second half + e = list(map(round, linspace(b/2, 0.0 , n-n2))) + for i in range(n2, n): + x[i] = uniform(-1.0, 1.0) * exp2(e[i - n2]) + y[i] = (uniform(-1.0, 1.0) * exp2(e[i - n2]) - DotExact(x, y)) / x[i] + + # Shuffle + pairs = list(zip(x, y)) + shuffle(pairs) + x, y = zip(*pairs) + + return DotExample(x, y, DotExact(x, y), Condition(x, y)) + + def RelativeError(res, ex): + x, y, target_sumprod, condition = ex + n = DotExact(list(x) + [-res], list(y) + [1]) + return fabs(n / target_sumprod) + + def Trial(dotfunc, c, n): + ex = GenDot(10, c) + res = dotfunc(ex.x, ex.y) + return RelativeError(res, ex) + + times = 1000 # Number of trials + n = 20 # Length of vectors + c = 1e30 # Target condition number + + # If the following test fails, it means that the C math library + # implementation of fma() is not compliant with the C99 standard + # and is inaccurate. To solve this problem, make a new build + # with the symbol UNRELIABLE_FMA defined. That will enable a + # slower but accurate code path that avoids the fma() call. + relative_err = median(Trial(math.sumprod, c, n) for i in range(times)) + self.assertLess(relative_err, 1e-16) + def testModf(self): self.assertRaises(TypeError, math.modf) @@ -1235,6 +1570,7 @@ def testPow(self): self.assertTrue(math.isnan(math.pow(2, NAN))) self.assertTrue(math.isnan(math.pow(0, NAN))) self.assertEqual(math.pow(1, NAN), 1) + self.assertRaises(OverflowError, math.pow, 1e+100, 1e+100) # pow(0., x) self.assertEqual(math.pow(0., INF), 0.) @@ -1550,7 +1886,7 @@ def testTan(self): try: self.assertTrue(math.isnan(math.tan(INF))) self.assertTrue(math.isnan(math.tan(NINF))) - except: + except ValueError: self.assertRaises(ValueError, math.tan, INF) self.assertRaises(ValueError, math.tan, NINF) self.assertTrue(math.isnan(math.tan(NAN))) @@ -1591,6 +1927,8 @@ def __trunc__(self): return 23 class TestNoTrunc: pass + class TestBadTrunc: + __trunc__ = BadDescr() self.assertEqual(math.trunc(TestTrunc()), 23) self.assertEqual(math.trunc(FloatTrunc()), 23) @@ -1599,6 +1937,7 @@ class TestNoTrunc: self.assertRaises(TypeError, math.trunc, 1, 2) self.assertRaises(TypeError, math.trunc, FloatLike(23.5)) self.assertRaises(TypeError, math.trunc, TestNoTrunc()) + self.assertRaises(ValueError, math.trunc, TestBadTrunc()) def testIsfinite(self): self.assertTrue(math.isfinite(0.0)) @@ -1626,11 +1965,11 @@ def testIsinf(self): self.assertFalse(math.isinf(0.)) self.assertFalse(math.isinf(1.)) - @requires_IEEE_754 def test_nan_constant(self): + # `math.nan` must be a quiet NaN with positive sign bit self.assertTrue(math.isnan(math.nan)) + self.assertEqual(math.copysign(1., math.nan), 1.) - @requires_IEEE_754 def test_inf_constant(self): self.assertTrue(math.isinf(math.inf)) self.assertGreater(math.inf, 0.0) @@ -1719,6 +2058,13 @@ def test_testfile(self): except OverflowError: result = 'OverflowError' + # C99+ says for math.h's sqrt: If the argument is +∞ or ±0, it is + # returned, unmodified. On another hand, for csqrt: If z is ±0+0i, + # the result is +0+0i. Lets correct zero sign of er to follow + # first convention. + if id in ['sqrt0002', 'sqrt0003', 'sqrt1001', 'sqrt1023']: + er = math.copysign(er, ar) + # Default tolerances ulp_tol, abs_tol = 5, 0.0 @@ -1802,6 +2148,8 @@ def test_mtestfile(self): '\n '.join(failures)) def test_prod(self): + from fractions import Fraction as F + prod = math.prod self.assertEqual(prod([]), 1) self.assertEqual(prod([], start=5), 5) @@ -1813,6 +2161,14 @@ def test_prod(self): self.assertEqual(prod([1.0, 2.0, 3.0, 4.0, 5.0]), 120.0) self.assertEqual(prod([1, 2, 3, 4.0, 5.0]), 120.0) self.assertEqual(prod([1.0, 2.0, 3.0, 4, 5]), 120.0) + self.assertEqual(prod([1., F(3, 2)]), 1.5) + + # Error in multiplication + class BadMultiply: + def __rmul__(self, other): + raise RuntimeError + with self.assertRaises(RuntimeError): + prod([10., BadMultiply()]) # Test overflow in fast-path for integers self.assertEqual(prod([1, 1, 2**32, 1, 1]), 2**32) @@ -2044,11 +2400,20 @@ def test_nextafter(self): float.fromhex('0x1.fffffffffffffp-1')) self.assertEqual(math.nextafter(1.0, INF), float.fromhex('0x1.0000000000001p+0')) + self.assertEqual(math.nextafter(1.0, -INF, steps=1), + float.fromhex('0x1.fffffffffffffp-1')) + self.assertEqual(math.nextafter(1.0, INF, steps=1), + float.fromhex('0x1.0000000000001p+0')) + self.assertEqual(math.nextafter(1.0, -INF, steps=3), + float.fromhex('0x1.ffffffffffffdp-1')) + self.assertEqual(math.nextafter(1.0, INF, steps=3), + float.fromhex('0x1.0000000000003p+0')) # x == y: y is returned - self.assertEqual(math.nextafter(2.0, 2.0), 2.0) - self.assertEqualSign(math.nextafter(-0.0, +0.0), +0.0) - self.assertEqualSign(math.nextafter(+0.0, -0.0), -0.0) + for steps in range(1, 5): + self.assertEqual(math.nextafter(2.0, 2.0, steps=steps), 2.0) + self.assertEqualSign(math.nextafter(-0.0, +0.0, steps=steps), +0.0) + self.assertEqualSign(math.nextafter(+0.0, -0.0, steps=steps), -0.0) # around 0.0 smallest_subnormal = sys.float_info.min * sys.float_info.epsilon @@ -2073,6 +2438,11 @@ def test_nextafter(self): self.assertIsNaN(math.nextafter(1.0, NAN)) self.assertIsNaN(math.nextafter(NAN, NAN)) + self.assertEqual(1.0, math.nextafter(1.0, INF, steps=0)) + with self.assertRaises(ValueError): + math.nextafter(1.0, INF, steps=-1) + + @requires_IEEE_754 def test_ulp(self): self.assertEqual(math.ulp(1.0), sys.float_info.epsilon) @@ -2112,6 +2482,14 @@ def __float__(self): # argument to a float. self.assertFalse(getattr(y, "converted", False)) + def test_input_exceptions(self): + self.assertRaises(TypeError, math.exp, "spam") + self.assertRaises(TypeError, math.erf, "spam") + self.assertRaises(TypeError, math.atan2, "spam", 1.0) + self.assertRaises(TypeError, math.atan2, 1.0, "spam") + self.assertRaises(TypeError, math.atan2, 1.0) + self.assertRaises(TypeError, math.atan2, 1.0, 2.0, 3.0) + # Custom assertions. def assertIsNaN(self, value): @@ -2252,7 +2630,7 @@ def test_fractions(self): def load_tests(loader, tests, pattern): from doctest import DocFileSuite - # tests.addTest(DocFileSuite("ieee754.txt")) + tests.addTest(DocFileSuite("ieee754.txt")) return tests if __name__ == '__main__': diff --git a/common/src/float_ops.rs b/common/src/float_ops.rs index 69ae8833a2..8e16e71c14 100644 --- a/common/src/float_ops.rs +++ b/common/src/float_ops.rs @@ -133,6 +133,69 @@ pub fn nextafter(x: f64, y: f64) -> f64 { } } +#[allow(clippy::float_cmp)] +pub fn nextafter_with_steps(x: f64, y: f64, steps: u64) -> f64 { + if x == y { + y + } else if x.is_nan() || y.is_nan() { + f64::NAN + } else if x >= f64::INFINITY { + f64::MAX + } else if x <= f64::NEG_INFINITY { + f64::MIN + } else if x == 0.0 { + f64::from_bits(1).copysign(y) + } else { + if steps == 0 { + return x; + } + + if x.is_nan() { + return x; + } + + if y.is_nan() { + return y; + } + + let sign_bit: u64 = 1 << 63; + + let mut ux = x.to_bits(); + let uy = y.to_bits(); + + let ax = ux & !sign_bit; + let ay = uy & !sign_bit; + + // If signs are different + if ((ux ^ uy) & sign_bit) != 0 { + return if ax + ay <= steps { + f64::from_bits(uy) + } else if ax < steps { + let result = (uy & sign_bit) | (steps - ax); + f64::from_bits(result) + } else { + ux -= steps; + f64::from_bits(ux) + }; + } + + // If signs are the same + if ax > ay { + if ax - ay >= steps { + ux -= steps; + f64::from_bits(ux) + } else { + f64::from_bits(uy) + } + } else if ay - ax >= steps { + ux += steps; + f64::from_bits(ux) + } else { + f64::from_bits(uy) + } + } +} + pub fn ulp(x: f64) -> f64 { if x.is_nan() { return x; diff --git a/stdlib/src/math.rs b/stdlib/src/math.rs index 30b088f714..fc4e740e40 100644 --- a/stdlib/src/math.rs +++ b/stdlib/src/math.rs @@ -629,7 +629,7 @@ mod math { #[pyfunction] fn fsum(seq: ArgIterable, vm: &VirtualMachine) -> PyResult { - let mut partials = vec![]; + let mut partials = Vec::with_capacity(32); let mut special_sum = 0.0; let mut inf_sum = 0.0; @@ -637,11 +637,11 @@ mod math { let mut x = *obj?; let xsave = x; - let mut j = 0; + let mut i = 0; // This inner loop applies `hi`/`lo` summation to each // partial so that the list of partial sums remains exact. - for i in 0..partials.len() { - let mut y: f64 = partials[i]; + for j in 0..partials.len() { + let mut y: f64 = partials[j]; if x.abs() < y.abs() { std::mem::swap(&mut x, &mut y); } @@ -650,33 +650,33 @@ mod math { let hi = x + y; let lo = y - (hi - x); if lo != 0.0 { - partials[j] = lo; - j += 1; + partials[i] = lo; + i += 1; } x = hi; } - if !x.is_finite() { - // a nonfinite x could arise either as - // a result of intermediate overflow, or - // as a result of a nan or inf in the - // summands - if xsave.is_finite() { - return Err(vm.new_overflow_error("intermediate overflow in fsum".to_owned())); - } - if xsave.is_infinite() { - inf_sum += xsave; + partials.truncate(i); + if x != 0.0 { + if !x.is_finite() { + // a nonfinite x could arise either as + // a result of intermediate overflow, or + // as a result of a nan or inf in the + // summands + if xsave.is_finite() { + return Err( + vm.new_overflow_error("intermediate overflow in fsum".to_owned()) + ); + } + if xsave.is_infinite() { + inf_sum += xsave; + } + special_sum += xsave; + // reset partials + partials.clear(); + } else { + partials.push(x); } - special_sum += xsave; - // reset partials - partials.clear(); - } - - if j >= partials.len() { - partials.push(x); - } else { - partials[j] = x; - partials.truncate(j + 1); } } if special_sum != 0.0 { @@ -831,9 +831,38 @@ mod math { (x.fract(), x.trunc()) } - #[pyfunction] - fn nextafter(x: ArgIntoFloat, y: ArgIntoFloat) -> f64 { - float_ops::nextafter(*x, *y) + #[derive(FromArgs)] + struct NextAfterArgs { + #[pyarg(positional)] + x: ArgIntoFloat, + #[pyarg(positional)] + y: ArgIntoFloat, + #[pyarg(named, optional)] + steps: OptionalArg, + } + + #[pyfunction] + fn nextafter(arg: NextAfterArgs, vm: &VirtualMachine) -> PyResult { + let steps: Option = arg + .steps + .map(|v| v.try_to_primitive(vm)) + .transpose()? + .into_option(); + match steps { + Some(steps) => { + if steps < 0 { + return Err( + vm.new_value_error("steps must be a non-negative integer".to_string()) + ); + } + Ok(float_ops::nextafter_with_steps( + *arg.x, + *arg.y, + steps as u64, + )) + } + None => Ok(float_ops::nextafter(*arg.x, *arg.y)), + } } #[pyfunction] @@ -917,10 +946,38 @@ mod math { // refer: https://github.com/python/cpython/blob/main/Modules/mathmodule.c#L3093-L3193 for obj in iter.iter(vm)? { let obj = obj?; + result = vm._mul(&result, &obj)?; + } + + Ok(result) + } - result = vm - ._mul(&result, &obj) - .map_err(|_| vm.new_type_error("math type error".to_owned()))?; + #[pyfunction] + fn sumprod( + p: ArgIterable, + q: ArgIterable, + vm: &VirtualMachine, + ) -> PyResult { + let mut p_iter = p.iter(vm)?; + let mut q_iter = q.iter(vm)?; + // We cannot just create a float because the iterator may contain + // anything as long as it supports __add__ and __mul__. + let mut result = vm.new_pyobj(0); + loop { + let m_p = p_iter.next(); + let m_q = q_iter.next(); + match (m_p, m_q) { + (Some(r_p), Some(r_q)) => { + let p = r_p?; + let q = r_q?; + let tmp = vm._mul(&p, &q)?; + result = vm._add(&result, &tmp)?; + } + (None, None) => break, + _ => { + return Err(vm.new_value_error("Inputs are not the same length".to_string())); + } + } } Ok(result) diff --git a/vm/src/builtins/float.rs b/vm/src/builtins/float.rs index 1cd041b7b9..27fdff12ee 100644 --- a/vm/src/builtins/float.rs +++ b/vm/src/builtins/float.rs @@ -124,8 +124,8 @@ fn inner_divmod(v1: f64, v2: f64, vm: &VirtualMachine) -> PyResult<(f64, f64)> { pub fn float_pow(v1: f64, v2: f64, vm: &VirtualMachine) -> PyResult { if v1.is_zero() && v2.is_sign_negative() { - let msg = format!("{v1} cannot be raised to a negative power"); - Err(vm.new_zero_division_error(msg)) + let msg = "0.0 cannot be raised to a negative power"; + Err(vm.new_zero_division_error(msg.to_owned())) } else if v1.is_sign_negative() && (v2.floor() - v2).abs() > f64::EPSILON { let v1 = Complex64::new(v1, 0.); let v2 = Complex64::new(v2, 0.);