Skip to content

bpo-29282: Add math.fma(): fused multiply-add function #17987

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

bpo-29282: Add math.fma(): fused multiply-add function #17987

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions Doc/library/math.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,22 @@ Number-theoretic and representation functions
:class:`~numbers.Integral` value.


.. function:: fma(x, y, z)

Fused multiply-add operation. Return ``(x * y) + z``, computed as though with
infinite precision and range followed by a single round to the ``float``
format. This operation often provides better accuracy than the direct
expression ``(x * y) + z``.

This function follows the specification of the fusedMultiplyAdd operation
described in the IEEE 754-2008 standard. The standard leaves one case
implementation-defined, namely the result of ``fma(0, inf, nan)``
and ``fma(inf, 0, nan)``. In these cases, ``math.fma`` returns a NaN,
and does not raise any exception.

.. versionadded:: 3.9


.. function:: fmod(x, y)

Return ``fmod(x, y)``, as defined by the platform C library. Note that the
Expand Down
7 changes: 7 additions & 0 deletions Doc/whatsnew/3.9.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,13 @@ Add :func:`math.ulp`: return the value of the least significant bit
of a float.
(Contributed by Victor Stinner in :issue:`39310`.)

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

nntplib
-------

Expand Down
233 changes: 233 additions & 0 deletions Lib/test/test_math.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2087,11 +2087,244 @@ def test_ulp(self):
self.assertEqual(math.ulp(-x), math.ulp(x))


class FMATests(unittest.TestCase):
""" Tests for math.fma. """

def test_fma_nan_results(self):
# Selected representative values.
values = [
-math.inf, -1e300, -2.3, -1e-300, -0.0,
0.0, 1e-300, 2.3, 1e300, math.inf, math.nan
]

# If any input is a NaN, the result should be a NaN, too.
for a, b in itertools.product(values, repeat=2):
self.assertIsNaN(math.fma(math.nan, a, b))
self.assertIsNaN(math.fma(a, math.nan, b))
self.assertIsNaN(math.fma(a, b, math.nan))

def test_fma_infinities(self):
# Cases involving infinite inputs or results.
positives = [1e-300, 2.3, 1e300, math.inf]
finites = [-1e300, -2.3, -1e-300, -0.0, 0.0, 1e-300, 2.3, 1e300]
non_nans = [-math.inf, -2.3, -0.0, 0.0, 2.3, math.inf]

# ValueError due to inf * 0 computation.
for c in non_nans:
for infinity in [math.inf, -math.inf]:
for zero in [0.0, -0.0]:
with self.assertRaises(ValueError):
math.fma(infinity, zero, c)
with self.assertRaises(ValueError):
math.fma(zero, infinity, c)

# ValueError when a*b and c both infinite of opposite signs.
for b in positives:
with self.assertRaises(ValueError):
math.fma(math.inf, b, -math.inf)
with self.assertRaises(ValueError):
math.fma(math.inf, -b, math.inf)
with self.assertRaises(ValueError):
math.fma(-math.inf, -b, -math.inf)
with self.assertRaises(ValueError):
math.fma(-math.inf, b, math.inf)
with self.assertRaises(ValueError):
math.fma(b, math.inf, -math.inf)
with self.assertRaises(ValueError):
math.fma(-b, math.inf, math.inf)
with self.assertRaises(ValueError):
math.fma(-b, -math.inf, -math.inf)
with self.assertRaises(ValueError):
math.fma(b, -math.inf, math.inf)

# Infinite result when a*b and c both infinite of the same sign.
for b in positives:
self.assertEqual(math.fma(math.inf, b, math.inf), math.inf)
self.assertEqual(math.fma(math.inf, -b, -math.inf), -math.inf)
self.assertEqual(math.fma(-math.inf, -b, math.inf), math.inf)
self.assertEqual(math.fma(-math.inf, b, -math.inf), -math.inf)
self.assertEqual(math.fma(b, math.inf, math.inf), math.inf)
self.assertEqual(math.fma(-b, math.inf, -math.inf), -math.inf)
self.assertEqual(math.fma(-b, -math.inf, math.inf), math.inf)
self.assertEqual(math.fma(b, -math.inf, -math.inf), -math.inf)

# Infinite result when a*b finite, c infinite.
for a, b in itertools.product(finites, finites):
self.assertEqual(math.fma(a, b, math.inf), math.inf)
self.assertEqual(math.fma(a, b, -math.inf), -math.inf)

# Infinite result when a*b infinite, c finite.
for b, c in itertools.product(positives, finites):
self.assertEqual(math.fma(math.inf, b, c), math.inf)
self.assertEqual(math.fma(-math.inf, b, c), -math.inf)
self.assertEqual(math.fma(-math.inf, -b, c), math.inf)
self.assertEqual(math.fma(math.inf, -b, c), -math.inf)

self.assertEqual(math.fma(b, math.inf, c), math.inf)
self.assertEqual(math.fma(b, -math.inf, c), -math.inf)
self.assertEqual(math.fma(-b, -math.inf, c), math.inf)
self.assertEqual(math.fma(-b, math.inf, c), -math.inf)

def test_fma_zero_result(self):
nonnegative_finites = [0.0, 1e-300, 2.3, 1e300]

# Zero results from exact zero inputs.
for b in nonnegative_finites:
self.assertIsPositiveZero(math.fma(0.0, b, 0.0))
self.assertIsPositiveZero(math.fma(0.0, b, -0.0))
self.assertIsNegativeZero(math.fma(0.0, -b, -0.0))
self.assertIsPositiveZero(math.fma(0.0, -b, 0.0))
self.assertIsPositiveZero(math.fma(-0.0, -b, 0.0))
self.assertIsPositiveZero(math.fma(-0.0, -b, -0.0))
self.assertIsNegativeZero(math.fma(-0.0, b, -0.0))
self.assertIsPositiveZero(math.fma(-0.0, b, 0.0))

self.assertIsPositiveZero(math.fma(b, 0.0, 0.0))
self.assertIsPositiveZero(math.fma(b, 0.0, -0.0))
self.assertIsNegativeZero(math.fma(-b, 0.0, -0.0))
self.assertIsPositiveZero(math.fma(-b, 0.0, 0.0))
self.assertIsPositiveZero(math.fma(-b, -0.0, 0.0))
self.assertIsPositiveZero(math.fma(-b, -0.0, -0.0))
self.assertIsNegativeZero(math.fma(b, -0.0, -0.0))
self.assertIsPositiveZero(math.fma(b, -0.0, 0.0))

# Exact zero result from nonzero inputs.
self.assertIsPositiveZero(math.fma(2.0, 2.0, -4.0))
self.assertIsPositiveZero(math.fma(2.0, -2.0, 4.0))
self.assertIsPositiveZero(math.fma(-2.0, -2.0, -4.0))
self.assertIsPositiveZero(math.fma(-2.0, 2.0, 4.0))

# Underflow to zero.
tiny = 1e-300
self.assertIsPositiveZero(math.fma(tiny, tiny, 0.0))
self.assertIsNegativeZero(math.fma(tiny, -tiny, 0.0))
self.assertIsPositiveZero(math.fma(-tiny, -tiny, 0.0))
self.assertIsNegativeZero(math.fma(-tiny, tiny, 0.0))
self.assertIsPositiveZero(math.fma(tiny, tiny, -0.0))
self.assertIsNegativeZero(math.fma(tiny, -tiny, -0.0))
self.assertIsPositiveZero(math.fma(-tiny, -tiny, -0.0))
self.assertIsNegativeZero(math.fma(-tiny, tiny, -0.0))

# Corner case where rounding the multiplication would
# give the wrong result.
x = float.fromhex('0x1p-500')
y = float.fromhex('0x1p-550')
z = float.fromhex('0x1p-1000')
self.assertIsNegativeZero(math.fma(x-y, x+y, -z))
self.assertIsPositiveZero(math.fma(y-x, x+y, z))
self.assertIsNegativeZero(math.fma(y-x, -(x+y), -z))
self.assertIsPositiveZero(math.fma(x-y, -(x+y), z))

def test_fma_overflow(self):
a = b = float.fromhex('0x1p512')
c = float.fromhex('0x1p1023')
# Overflow from multiplication.
with self.assertRaises(OverflowError):
math.fma(a, b, 0.0)
self.assertEqual(math.fma(a, b/2.0, 0.0), c)
# Overflow from the addition.
with self.assertRaises(OverflowError):
math.fma(a, b/2.0, c)
# No overflow, even though a*b overflows a float.
self.assertEqual(math.fma(a, b, -c), c)

# Extreme case: a * b is exactly at the overflow boundary, so the
# tiniest offset makes a difference between overflow and a finite
# result.
a = float.fromhex('0x1.ffffffc000000p+511')
b = float.fromhex('0x1.0000002000000p+512')
c = float.fromhex('0x0.0000000000001p-1022')
with self.assertRaises(OverflowError):
math.fma(a, b, 0.0)
with self.assertRaises(OverflowError):
math.fma(a, b, c)
self.assertEqual(math.fma(a, b, -c),
float.fromhex('0x1.fffffffffffffp+1023'))

# Another extreme case: here a*b is about as large as possible subject
# to math.fma(a, b, c) being finite.
a = float.fromhex('0x1.ae565943785f9p+512')
b = float.fromhex('0x1.3094665de9db8p+512')
c = float.fromhex('0x1.fffffffffffffp+1023')
self.assertEqual(math.fma(a, b, -c), c)

def test_fma_single_round(self):
a = float.fromhex('0x1p-50')
self.assertEqual(math.fma(a - 1.0, a + 1.0, 1.0), a*a)

def test_random(self):
# A collection of randomly generated inputs for which the naive FMA
# (with two rounds) gives a different result from a singly-rounded FMA.

# tuples (a, b, c, expected)
test_values = [
('0x1.694adde428b44p-1', '0x1.371b0d64caed7p-1',
'0x1.f347e7b8deab8p-4', '0x1.19f10da56c8adp-1'),
('0x1.605401ccc6ad6p-2', '0x1.ce3a40bf56640p-2',
'0x1.96e3bf7bf2e20p-2', '0x1.1af6d8aa83101p-1'),
('0x1.e5abd653a67d4p-2', '0x1.a2e400209b3e6p-1',
'0x1.a90051422ce13p-1', '0x1.37d68cc8c0fbbp+0'),
('0x1.f94e8efd54700p-2', '0x1.123065c812cebp-1',
'0x1.458f86fb6ccd0p-1', '0x1.ccdcee26a3ff3p-1'),
('0x1.bd926f1eedc96p-1', '0x1.eee9ca68c5740p-1',
'0x1.960c703eb3298p-2', '0x1.3cdcfb4fdb007p+0'),
('0x1.27348350fbccdp-1', '0x1.3b073914a53f1p-1',
'0x1.e300da5c2b4cbp-1', '0x1.4c51e9a3c4e29p+0'),
('0x1.2774f00b3497bp-1', '0x1.7038ec336bff0p-2',
'0x1.2f6f2ccc3576bp-1', '0x1.99ad9f9c2688bp-1'),
('0x1.51d5a99300e5cp-1', '0x1.5cd74abd445a1p-1',
'0x1.8880ab0bbe530p-1', '0x1.3756f96b91129p+0'),
('0x1.73cb965b821b8p-2', '0x1.218fd3d8d5371p-1',
'0x1.d1ea966a1f758p-2', '0x1.5217b8fd90119p-1'),
('0x1.4aa98e890b046p-1', '0x1.954d85dff1041p-1',
'0x1.122b59317ebdfp-1', '0x1.0bf644b340cc5p+0'),
('0x1.e28f29e44750fp-1', '0x1.4bcc4fdcd18fep-1',
'0x1.fd47f81298259p-1', '0x1.9b000afbc9995p+0'),
('0x1.d2e850717fe78p-3', '0x1.1dd7531c303afp-1',
'0x1.e0869746a2fc2p-2', '0x1.316df6eb26439p-1'),
('0x1.cf89c75ee6fbap-2', '0x1.b23decdc66825p-1',
'0x1.3d1fe76ac6168p-1', '0x1.00d8ea4c12abbp+0'),
('0x1.3265ae6f05572p-2', '0x1.16d7ec285f7a2p-1',
'0x1.0b8405b3827fbp-1', '0x1.5ef33c118a001p-1'),
('0x1.c4d1bf55ec1a5p-1', '0x1.bc59618459e12p-2',
'0x1.ce5b73dc1773dp-1', '0x1.496cf6164f99bp+0'),
('0x1.d350026ac3946p-1', '0x1.9a234e149a68cp-2',
'0x1.f5467b1911fd6p-2', '0x1.b5cee3225caa5p-1'),
]
for a_hex, b_hex, c_hex, expected_hex in test_values:
a = float.fromhex(a_hex)
b = float.fromhex(b_hex)
c = float.fromhex(c_hex)
expected = float.fromhex(expected_hex)
self.assertEqual(math.fma(a, b, c), expected)
self.assertEqual(math.fma(b, a, c), expected)

# Custom assertions.
def assertIsNaN(self, value):
self.assertTrue(
math.isnan(value),
msg="Expected a NaN, got {!r}".format(value)
)

def assertIsPositiveZero(self, value):
self.assertTrue(
value == 0 and math.copysign(1, value) > 0,
msg="Expected a positive zero, got {!r}".format(value)
)

def assertIsNegativeZero(self, value):
self.assertTrue(
value == 0 and math.copysign(1, value) < 0,
msg="Expected a negative zero, got {!r}".format(value)
)


def test_main():
from doctest import DocFileSuite
suite = unittest.TestSuite()
suite.addTest(unittest.makeSuite(MathTests))
suite.addTest(unittest.makeSuite(IsCloseTests))
suite.addTest(unittest.makeSuite(FMATests))
suite.addTest(DocFileSuite("ieee754.txt"))
run_unittest(suite)

Expand Down
2 changes: 2 additions & 0 deletions Misc/NEWS.d/next/Library/2020-01-13-15-04-26.bpo-29282.lO827U.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Added new :func:`math.fma` function, wrapping C99's ``fma()`` operation: fused
multiply-add function.
61 changes: 60 additions & 1 deletion Modules/clinic/mathmodule.c.h
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading