This PR is a copy the following https://hg.python.org/cpython/rev/b33012ef1417 written by Mark Dickinson. This change has been reverted, see: https://bugs.python.org/issue29282#msg285956

I created this PR to reconsider math.fma() and to use this PR as a starting point to fix the implementation. I'm now waiting the CIs to see how things are going on.

If tests continue to fail on some platforms, I plan to manually handle NaN and INF in the C code, before calling libc fma().

I would prefer to get the same behavior for fma(x, y, z) than (x * y) + z if possible. Here I don't understand why math.fma(math.inf, 0, 0) doesn't return nan. Is it a requirement of IEEE 754, or a Python deliberate choice?

>>> math.fma(math.inf, 0, math.nan)
nan
>>> math.fma(math.inf, 0, 0)
ValueError: invalid operation in fma

>>> (math.inf * 0) + math.nan
nan
>>> (math.inf * 0) + 0
nan

cc @mdickinson

🤖 New build scheduled with the buildbot fleet by @vstinner for commit 17d570c 🤖

If you want to schedule another build, you need to add the "🔨 test-with-buildbots" label again.

Strange. Tests passed on Linux and Windows pre-commit CIs. Let me try on stable buildbots.

Is it a requirement of IEEE 754, or a Python deliberate choice?

fma is a standard IEEE 754 function, so yes, the behaviour in corner cases is fully specified by IEEE 754, and we should expect to follow that behaviour.

Here I don't understand why math.fma(math.inf, 0, 0) doesn't return nan.

This is the usual rule about propagating NaNs, together with Python's rules for wrapping existing functions (as articulated at the top of mathmodule.c). In situations where IEEE 754 signals a floating-point exception, Python should be raising. But if any of the three arguments to fma is a NaN, then a NaN is silently propagated, and no floating-point exception is signaled.

Think of fma as a 3-argument operation, not as exactly equivalent to a*b + c.

Think of fma as a 3-argument operation, not as exactly equivalent to a*b + c.

You can also look at Decimal's implementation, which should also be following IEEE 754 rules (but with the added complication of having to handle signalling NaNs correctly).

test_fma_zero_result() fails on FreeBSD:

• AMD64 FreeBSD Non-Debug PR: https://buildbot.python.org/all/#/builders/291/builds/12
• AMD64 FreeBSD Shared PR: https://buildbot.python.org/all/#/builders/287/builds/12

FAIL: test_fma_zero_result (test.test_math.FMATests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/usr/home/buildbot/python/pull_request.koobs-freebsd-9e36.nondebug/build/Lib/test/test_math.py", line 2213, in test_fma_zero_result
    self.assertIsNegativeZero(math.fma(x-y, x+y, -z))
  File "/usr/home/buildbot/python/pull_request.koobs-freebsd-9e36.nondebug/build/Lib/test/test_math.py", line 2316, in assertIsNegativeZero
    self.assertTrue(
AssertionError: False is not true : Expected a negative zero, got 0.0

Shared pythoninfo:

CC.version: FreeBSD clang version 9.0.1 (git@github.com:llvm/llvm-project.git c1a0a213378a458fbea1a5c77b315c7dce08fd05) (based on LLVM 9.0.1)
...
platform.libc_ver: libc 7
platform.platform: FreeBSD-13.0-CURRENT-amd64-64bit-ELF

Non-debug pythoninfo:

CC.version: FreeBSD clang version 8.0.1 (tags/RELEASE_801/final 366581) (based on LLVM 8.0.1)
...
platform.libc_ver: libc 7
platform.platform: FreeBSD-12.1-RELEASE-p1-amd64-64bit-ELF

Note: PPC64 Fedora PR failure is unrelated (test_distutils: https://bugs.python.org/issue39248).

Failures on x86 Windows7 PR (32-bit):

======================================================================
ERROR: test_fma_overflow (test.test_math.FMATests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "D:\cygwin\home\db3l\buildarea\pull_request.bolen-windows7\build\lib\test\test_math.py", line 2241, in test_fma_overflow
    self.assertEqual(math.fma(a, b, -c),
OverflowError: overflow in fma

======================================================================
FAIL: test_random (test.test_math.FMATests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "D:\cygwin\home\db3l\buildarea\pull_request.bolen-windows7\build\lib\test\test_math.py", line 2299, in test_random
    self.assertEqual(math.fma(a, b, c), expected)
AssertionError: 0.5506672157701096 != 0.5506672157701097

pythoninfo:

sys.version: 3.9.0a2+ (heads/refs/pull/17937/merge:e4f904e, Jan 13 2020, 10:46:12) [MSC v.1913 32 bit (Intel)]
sys.version_info: sys.version_info(major=3, minor=9, micro=0, releaselevel='alpha', serial=2)
sys.windowsversion: sys.getwindowsversion(major=6, minor=1, build=7601, platform=2, service_pack='Service Pack 1')
sys.winver: 3.9-32

"MSC v.1913" is Visual Studio 2017 Update 6.

fma is a standard IEEE 754 function, so yes, the behaviour in corner cases is fully specified by IEEE 754

Sorry, I'm wrong: "fully" is incorrect. There's one case where both IEEE 754 and C99 refuse to specify the behaviour, namely fma(±0, ±inf, nan) and fma(±inf, ±0, nan). IEEE 754 says for these cases (see section 7.2(c) of IEEE 754-2019) that:

it is implementation defined whether the invalid operation exception is signaled

while C99 says in Annex F, 9.10.1:

fma(x, y, z) returns a NaN and optionally raises the invalid floating-point exception if one of x and y is infinite, the other is zero, and z is a NaN

(the emphasis above is mine). So for us, fma(0, inf, nan) could either raise ValueError or return a NaN. I don't know whether we should make a choice for Python in this case, or leave this system-defined.

My inclination would be to specify this in Python, and have Python return a NaN in this case, and not raise. @tim-one thoughts?

Failures on x86 Windows7 PR (32-bit):

Okay, so we still have the same issue as last time I tried this. :-(

I'm really not comfortable with knowingly delivering a substandard fma on a mainstream platform.

So for us, fma(0, inf, nan) could either raise ValueError or return a NaN. ...

My inclination would be to specify this in Python, and have Python return a NaN in this case, and not raise. @tim-one thoughts?

Yes, just return a NaN. anything + NaN is going to return a NaN, so complaining about what happens during the computation of anything is pointless.

About the failure of single-rounding on Windows. it's surprising, and a few minutes on Google only found complaints about this in the BPO issue report. I wonder whether we're using (or failing to use) some goofy Windows-specific compiler (linker?) flag(s), but can't make time now to thrash with that.

Failures on x86 Windows7 PR (32-bit):

Is that the only Windows build with a failure? Windows 7 end-of-life is tomorrow (14 Jan 2020), and nobody cares about 32-bit boxes anymore anyway 😉.

I checked out this branch and tried it on my own Windows box (Win 10 Pro, Visual Studo 2017). test_math passed under all 4 combinations of {Release, Debug} x {64-bit, 32-bit}.

Here's a sample ID line:

Python 3.9.0a2+ (heads/fma:17d570cb0c, Jan 13 2020, 19:34:03) [MSC v.1916 32 bit (Intel)] on win32

So the Win32 buildbot failure may be spurious. No idea why. Ancient OS? cygwin mucking with something? Slightly earlier build of Visual Studio? The buildbot is actually running on a 32-bit box? ...

Ah! I see now that the other Windows buildbots passed (as did my own Windows box), including an AMD64 Windows 7 SP1 bot.

So there's "something wrong" with the sole Windows oddball that failed. But I don't know how to pursue that.

Windows 7 end-of-life is tomorrow (14 Jan 2020), and nobody cares about 32-bit boxes anymore anyway 😉.

I wish that were true, but interactions with customers say otherwise. :-(

OTOH, none of those customers are going to be using Python 3.9 any time soon (one of them "upgraded" last year from Python 2 to Python 3.4).

So there's "something wrong" with the sole Windows oddball that failed.

My suspicion is that use of the x87 FPU plays some role in this.