Skip to content

Commit 17d570c

Browse files
vstinnermdickinson
vstinner
and
mdickinson
committed
bpo-29282: Add math.fma(): fused multiply-add function
Co-Authored-By: Mark Dickinson <mdickinson@enthought.com>
1 parent e7c9f4a commit 17d570c

File tree

6 files changed

+360
-1
lines changed

6 files changed

+360
-1
lines changed

‎Doc/library/math.rst

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,22 @@ Number-theoretic and representation functions
8282
:class:`~numbers.Integral` value.
8383

8484

85+
.. function:: fma(x, y, z)
86+
87+
Fused multiply-add operation. Return ``(x * y) + z``, computed as though with
88+
infinite precision and range followed by a single round to the ``float``
89+
format. This operation often provides better accuracy than the direct
90+
expression ``(x * y) + z``.
91+
92+
This function follows the specification of the fusedMultiplyAdd operation
93+
described in the IEEE 754-2008 standard. The standard leaves one case
94+
implementation-defined, namely the result of ``fma(0, inf, nan)``
95+
and ``fma(inf, 0, nan)``. In these cases, ``math.fma`` returns a NaN,
96+
and does not raise any exception.
97+
98+
.. versionadded:: 3.9
99+
100+
85101
.. function:: fmod(x, y)
86102

87103
Return ``fmod(x, y)``, as defined by the platform C library. Note that the

‎Doc/whatsnew/3.9.rst

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -188,6 +188,13 @@ Add :func:`math.ulp`: return the value of the least significant bit
188188
of a float.
189189
(Contributed by Victor Stinner in :issue:`39310`.)
190190

191+
A new function :func:`~math.fma` for fused multiply-add operations has been
192+
added. This function computes ``x * y + z`` with only a single round, and so
193+
avoids any intermediate loss of precision. It wraps the ``fma`` function
194+
provided by C99, and follows the specification of the IEEE 754-2008
195+
"fusedMultiplyAdd" operation for special cases.
196+
(Contributed by Mark Dickinson in :issue:`29282`.)
197+
191198
nntplib
192199
-------
193200

‎Lib/test/test_math.py

Lines changed: 233 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2087,11 +2087,244 @@ def test_ulp(self):
20872087
self.assertEqual(math.ulp(-x), math.ulp(x))
20882088

20892089

2090+
class FMATests(unittest.TestCase):
2091+
""" Tests for math.fma. """
2092+
2093+
def test_fma_nan_results(self):
2094+
# Selected representative values.
2095+
values = [
2096+
-math.inf, -1e300, -2.3, -1e-300, -0.0,
2097+
0.0, 1e-300, 2.3, 1e300, math.inf, math.nan
2098+
]
2099+
2100+
# If any input is a NaN, the result should be a NaN, too.
2101+
for a, b in itertools.product(values, repeat=2):
2102+
self.assertIsNaN(math.fma(math.nan, a, b))
2103+
self.assertIsNaN(math.fma(a, math.nan, b))
2104+
self.assertIsNaN(math.fma(a, b, math.nan))
2105+
2106+
def test_fma_infinities(self):
2107+
# Cases involving infinite inputs or results.
2108+
positives = [1e-300, 2.3, 1e300, math.inf]
2109+
finites = [-1e300, -2.3, -1e-300, -0.0, 0.0, 1e-300, 2.3, 1e300]
2110+
non_nans = [-math.inf, -2.3, -0.0, 0.0, 2.3, math.inf]
2111+
2112+
# ValueError due to inf * 0 computation.
2113+
for c in non_nans:
2114+
for infinity in [math.inf, -math.inf]:
2115+
for zero in [0.0, -0.0]:
2116+
with self.assertRaises(ValueError):
2117+
math.fma(infinity, zero, c)
2118+
with self.assertRaises(ValueError):
2119+
math.fma(zero, infinity, c)
2120+
2121+
# ValueError when a*b and c both infinite of opposite signs.
2122+
for b in positives:
2123+
with self.assertRaises(ValueError):
2124+
math.fma(math.inf, b, -math.inf)
2125+
with self.assertRaises(ValueError):
2126+
math.fma(math.inf, -b, math.inf)
2127+
with self.assertRaises(ValueError):
2128+
math.fma(-math.inf, -b, -math.inf)
2129+
with self.assertRaises(ValueError):
2130+
math.fma(-math.inf, b, math.inf)
2131+
with self.assertRaises(ValueError):
2132+
math.fma(b, math.inf, -math.inf)
2133+
with self.assertRaises(ValueError):
2134+
math.fma(-b, math.inf, math.inf)
2135+
with self.assertRaises(ValueError):
2136+
math.fma(-b, -math.inf, -math.inf)
2137+
with self.assertRaises(ValueError):
2138+
math.fma(b, -math.inf, math.inf)
2139+
2140+
# Infinite result when a*b and c both infinite of the same sign.
2141+
for b in positives:
2142+
self.assertEqual(math.fma(math.inf, b, math.inf), math.inf)
2143+
self.assertEqual(math.fma(math.inf, -b, -math.inf), -math.inf)
2144+
self.assertEqual(math.fma(-math.inf, -b, math.inf), math.inf)
2145+
self.assertEqual(math.fma(-math.inf, b, -math.inf), -math.inf)
2146+
self.assertEqual(math.fma(b, math.inf, math.inf), math.inf)
2147+
self.assertEqual(math.fma(-b, math.inf, -math.inf), -math.inf)
2148+
self.assertEqual(math.fma(-b, -math.inf, math.inf), math.inf)
2149+
self.assertEqual(math.fma(b, -math.inf, -math.inf), -math.inf)
2150+
2151+
# Infinite result when a*b finite, c infinite.
2152+
for a, b in itertools.product(finites, finites):
2153+
self.assertEqual(math.fma(a, b, math.inf), math.inf)
2154+
self.assertEqual(math.fma(a, b, -math.inf), -math.inf)
2155+
2156+
# Infinite result when a*b infinite, c finite.
2157+
for b, c in itertools.product(positives, finites):
2158+
self.assertEqual(math.fma(math.inf, b, c), math.inf)
2159+
self.assertEqual(math.fma(-math.inf, b, c), -math.inf)
2160+
self.assertEqual(math.fma(-math.inf, -b, c), math.inf)
2161+
self.assertEqual(math.fma(math.inf, -b, c), -math.inf)
2162+
2163+
self.assertEqual(math.fma(b, math.inf, c), math.inf)
2164+
self.assertEqual(math.fma(b, -math.inf, c), -math.inf)
2165+
self.assertEqual(math.fma(-b, -math.inf, c), math.inf)
2166+
self.assertEqual(math.fma(-b, math.inf, c), -math.inf)
2167+
2168+
def test_fma_zero_result(self):
2169+
nonnegative_finites = [0.0, 1e-300, 2.3, 1e300]
2170+
2171+
# Zero results from exact zero inputs.
2172+
for b in nonnegative_finites:
2173+
self.assertIsPositiveZero(math.fma(0.0, b, 0.0))
2174+
self.assertIsPositiveZero(math.fma(0.0, b, -0.0))
2175+
self.assertIsNegativeZero(math.fma(0.0, -b, -0.0))
2176+
self.assertIsPositiveZero(math.fma(0.0, -b, 0.0))
2177+
self.assertIsPositiveZero(math.fma(-0.0, -b, 0.0))
2178+
self.assertIsPositiveZero(math.fma(-0.0, -b, -0.0))
2179+
self.assertIsNegativeZero(math.fma(-0.0, b, -0.0))
2180+
self.assertIsPositiveZero(math.fma(-0.0, b, 0.0))
2181+
2182+
self.assertIsPositiveZero(math.fma(b, 0.0, 0.0))
2183+
self.assertIsPositiveZero(math.fma(b, 0.0, -0.0))
2184+
self.assertIsNegativeZero(math.fma(-b, 0.0, -0.0))
2185+
self.assertIsPositiveZero(math.fma(-b, 0.0, 0.0))
2186+
self.assertIsPositiveZero(math.fma(-b, -0.0, 0.0))
2187+
self.assertIsPositiveZero(math.fma(-b, -0.0, -0.0))
2188+
self.assertIsNegativeZero(math.fma(b, -0.0, -0.0))
2189+
self.assertIsPositiveZero(math.fma(b, -0.0, 0.0))
2190+
2191+
# Exact zero result from nonzero inputs.
2192+
self.assertIsPositiveZero(math.fma(2.0, 2.0, -4.0))
2193+
self.assertIsPositiveZero(math.fma(2.0, -2.0, 4.0))
2194+
self.assertIsPositiveZero(math.fma(-2.0, -2.0, -4.0))
2195+
self.assertIsPositiveZero(math.fma(-2.0, 2.0, 4.0))
2196+
2197+
# Underflow to zero.
2198+
tiny = 1e-300
2199+
self.assertIsPositiveZero(math.fma(tiny, tiny, 0.0))
2200+
self.assertIsNegativeZero(math.fma(tiny, -tiny, 0.0))
2201+
self.assertIsPositiveZero(math.fma(-tiny, -tiny, 0.0))
2202+
self.assertIsNegativeZero(math.fma(-tiny, tiny, 0.0))
2203+
self.assertIsPositiveZero(math.fma(tiny, tiny, -0.0))
2204+
self.assertIsNegativeZero(math.fma(tiny, -tiny, -0.0))
2205+
self.assertIsPositiveZero(math.fma(-tiny, -tiny, -0.0))
2206+
self.assertIsNegativeZero(math.fma(-tiny, tiny, -0.0))
2207+
2208+
# Corner case where rounding the multiplication would
2209+
# give the wrong result.
2210+
x = float.fromhex('0x1p-500')
2211+
y = float.fromhex('0x1p-550')
2212+
z = float.fromhex('0x1p-1000')
2213+
self.assertIsNegativeZero(math.fma(x-y, x+y, -z))
2214+
self.assertIsPositiveZero(math.fma(y-x, x+y, z))
2215+
self.assertIsNegativeZero(math.fma(y-x, -(x+y), -z))
2216+
self.assertIsPositiveZero(math.fma(x-y, -(x+y), z))
2217+
2218+
def test_fma_overflow(self):
2219+
a = b = float.fromhex('0x1p512')
2220+
c = float.fromhex('0x1p1023')
2221+
# Overflow from multiplication.
2222+
with self.assertRaises(OverflowError):
2223+
math.fma(a, b, 0.0)
2224+
self.assertEqual(math.fma(a, b/2.0, 0.0), c)
2225+
# Overflow from the addition.
2226+
with self.assertRaises(OverflowError):
2227+
math.fma(a, b/2.0, c)
2228+
# No overflow, even though a*b overflows a float.
2229+
self.assertEqual(math.fma(a, b, -c), c)
2230+
2231+
# Extreme case: a * b is exactly at the overflow boundary, so the
2232+
# tiniest offset makes a difference between overflow and a finite
2233+
# result.
2234+
a = float.fromhex('0x1.ffffffc000000p+511')
2235+
b = float.fromhex('0x1.0000002000000p+512')
2236+
c = float.fromhex('0x0.0000000000001p-1022')
2237+
with self.assertRaises(OverflowError):
2238+
math.fma(a, b, 0.0)
2239+
with self.assertRaises(OverflowError):
2240+
math.fma(a, b, c)
2241+
self.assertEqual(math.fma(a, b, -c),
2242+
float.fromhex('0x1.fffffffffffffp+1023'))
2243+
2244+
# Another extreme case: here a*b is about as large as possible subject
2245+
# to math.fma(a, b, c) being finite.
2246+
a = float.fromhex('0x1.ae565943785f9p+512')
2247+
b = float.fromhex('0x1.3094665de9db8p+512')
2248+
c = float.fromhex('0x1.fffffffffffffp+1023')
2249+
self.assertEqual(math.fma(a, b, -c), c)
2250+
2251+
def test_fma_single_round(self):
2252+
a = float.fromhex('0x1p-50')
2253+
self.assertEqual(math.fma(a - 1.0, a + 1.0, 1.0), a*a)
2254+
2255+
def test_random(self):
2256+
# A collection of randomly generated inputs for which the naive FMA
2257+
# (with two rounds) gives a different result from a singly-rounded FMA.
2258+
2259+
# tuples (a, b, c, expected)
2260+
test_values = [
2261+
('0x1.694adde428b44p-1', '0x1.371b0d64caed7p-1',
2262+
'0x1.f347e7b8deab8p-4', '0x1.19f10da56c8adp-1'),
2263+
('0x1.605401ccc6ad6p-2', '0x1.ce3a40bf56640p-2',
2264+
'0x1.96e3bf7bf2e20p-2', '0x1.1af6d8aa83101p-1'),
2265+
('0x1.e5abd653a67d4p-2', '0x1.a2e400209b3e6p-1',
2266+
'0x1.a90051422ce13p-1', '0x1.37d68cc8c0fbbp+0'),
2267+
('0x1.f94e8efd54700p-2', '0x1.123065c812cebp-1',
2268+
'0x1.458f86fb6ccd0p-1', '0x1.ccdcee26a3ff3p-1'),
2269+
('0x1.bd926f1eedc96p-1', '0x1.eee9ca68c5740p-1',
2270+
'0x1.960c703eb3298p-2', '0x1.3cdcfb4fdb007p+0'),
2271+
('0x1.27348350fbccdp-1', '0x1.3b073914a53f1p-1',
2272+
'0x1.e300da5c2b4cbp-1', '0x1.4c51e9a3c4e29p+0'),
2273+
('0x1.2774f00b3497bp-1', '0x1.7038ec336bff0p-2',
2274+
'0x1.2f6f2ccc3576bp-1', '0x1.99ad9f9c2688bp-1'),
2275+
('0x1.51d5a99300e5cp-1', '0x1.5cd74abd445a1p-1',
2276+
'0x1.8880ab0bbe530p-1', '0x1.3756f96b91129p+0'),
2277+
('0x1.73cb965b821b8p-2', '0x1.218fd3d8d5371p-1',
2278+
'0x1.d1ea966a1f758p-2', '0x1.5217b8fd90119p-1'),
2279+
('0x1.4aa98e890b046p-1', '0x1.954d85dff1041p-1',
2280+
'0x1.122b59317ebdfp-1', '0x1.0bf644b340cc5p+0'),
2281+
('0x1.e28f29e44750fp-1', '0x1.4bcc4fdcd18fep-1',
2282+
'0x1.fd47f81298259p-1', '0x1.9b000afbc9995p+0'),
2283+
('0x1.d2e850717fe78p-3', '0x1.1dd7531c303afp-1',
2284+
'0x1.e0869746a2fc2p-2', '0x1.316df6eb26439p-1'),
2285+
('0x1.cf89c75ee6fbap-2', '0x1.b23decdc66825p-1',
2286+
'0x1.3d1fe76ac6168p-1', '0x1.00d8ea4c12abbp+0'),
2287+
('0x1.3265ae6f05572p-2', '0x1.16d7ec285f7a2p-1',
2288+
'0x1.0b8405b3827fbp-1', '0x1.5ef33c118a001p-1'),
2289+
('0x1.c4d1bf55ec1a5p-1', '0x1.bc59618459e12p-2',
2290+
'0x1.ce5b73dc1773dp-1', '0x1.496cf6164f99bp+0'),
2291+
('0x1.d350026ac3946p-1', '0x1.9a234e149a68cp-2',
2292+
'0x1.f5467b1911fd6p-2', '0x1.b5cee3225caa5p-1'),
2293+
]
2294+
for a_hex, b_hex, c_hex, expected_hex in test_values:
2295+
a = float.fromhex(a_hex)
2296+
b = float.fromhex(b_hex)
2297+
c = float.fromhex(c_hex)
2298+
expected = float.fromhex(expected_hex)
2299+
self.assertEqual(math.fma(a, b, c), expected)
2300+
self.assertEqual(math.fma(b, a, c), expected)
2301+
2302+
# Custom assertions.
2303+
def assertIsNaN(self, value):
2304+
self.assertTrue(
2305+
math.isnan(value),
2306+
msg="Expected a NaN, got {!r}".format(value)
2307+
)
2308+
2309+
def assertIsPositiveZero(self, value):
2310+
self.assertTrue(
2311+
value == 0 and math.copysign(1, value) > 0,
2312+
msg="Expected a positive zero, got {!r}".format(value)
2313+
)
2314+
2315+
def assertIsNegativeZero(self, value):
2316+
self.assertTrue(
2317+
value == 0 and math.copysign(1, value) < 0,
2318+
msg="Expected a negative zero, got {!r}".format(value)
2319+
)
2320+
2321+
20902322
def test_main():
20912323
from doctest import DocFileSuite
20922324
suite = unittest.TestSuite()
20932325
suite.addTest(unittest.makeSuite(MathTests))
20942326
suite.addTest(unittest.makeSuite(IsCloseTests))
2327+
suite.addTest(unittest.makeSuite(FMATests))
20952328
suite.addTest(DocFileSuite("ieee754.txt"))
20962329
run_unittest(suite)
20972330

‎Misc/NEWS.d/next/Library/2020-01-13-15-04-26.bpo-29282.lO827U.rst

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
Added new :func:`math.fma` function, wrapping C99's ``fma()`` operation: fused
2+
multiply-add function.

‎Modules/clinic/mathmodule.c.h

Lines changed: 60 additions & 1 deletion
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)