Fused multiply-add: proposal to add math.fma()

@mdickinson

BPO	29282
Nosy	@mdickinson, @pitrou, @vstinner, @stevendaprano, @skrah, @nschloe, @serhiy-storchaka, @efahl, @cdonovick
PRs	bpo-29282: Add math.fma(): fused multiply-add function #17987
Files	xmathmodule.c example.py setup.py fma_reference.py math_fma.patch math_fma2.patch math_fma3.patch

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = None
closed_at = None
created_at = <Date 2017-01-16.10:01:37.166>
labels = ['3.7', 'type-feature', 'library']
title = 'Fused multiply-add: proposal to add math.fma()'
updated_at = <Date 2021-06-19.09:24:45.622>
user = 'https://bugs.python.org/jurajsukop'

bugs.python.org fields:

activity = <Date 2021-06-19.09:24:45.622>
actor = 'mark.dickinson'
assignee = 'none'
closed = False
closed_date = None
closer = None
components = ['Library (Lib)']
creation = <Date 2017-01-16.10:01:37.166>
creator = 'juraj.sukop'
dependencies = []
files = ['46298', '46299', '46300', '46304', '46321', '46345', '46346']
hgrepos = []
issue_num = 29282
keywords = ['patch']
message_count = 36.0
messages = ['285547', '285549', '285550', '285551', '285552', '285553', '285557', '285564', '285565', '285566', '285582', '285682', '285685', '285836', '285837', '285839', '285840', '285841', '285843', '285844', '285955', '285956', '285957', '285958', '285959', '286201', '312431', '312432', '312433', '312480', '312488', '359903', '359911', '359913', '359921', '359926']
nosy_count = 11.0
nosy_names = ['mark.dickinson', 'pitrou', 'vstinner', 'steven.daprano', 'skrah', 'nschloe', 'python-dev', 'serhiy.storchaka', 'eric.fahlgren', 'juraj.sukop', 'donovick']
pr_nums = ['17987']
priority = 'normal'
resolution = None
stage = 'patch review'
status = 'open'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue29282'
versions = ['Python 3.7']

Linked PRs

gh-73468: Add math.fma() function #116667

Fused multiply-add (henceforth FMA) is an operation which calculates the product of two numbers and then the sum of the product and a third number with just one floating-point rounding. More concretely:

    r = x*y + z

The value of r is the same as if the RHS was calculated with infinite precision and the rounded to a 32-bit single-precision or 64-bit double-precision floating-point number [1].

Even though one FMA CPU instruction might be calculated faster than the two separate instructions for multiply and add, its main advantage comes from the increased precision of numerical computations that involve the accumulation of products. Examples which benefit from using FMA are: dot product [2], compensated arithmetic [3], polynomial evaluation [4], matrix multiplication, Newton's method and many more [5].

C99 includes [6] fma function to math.h and emulates the calculation if the FMA instruction is not present on the host CPU [7]. PEP-7 states that "Python versions greater than or equal to 3.6 use C89 with several select C99 features" and that "Future C99 features may be added to this list in the future depending on compiler support" [8].

This proposal is then about adding new fma function with the following signature to math module:

    math.fma(x, y, z)
    
'''Return a float representing the result of the operation `x*y + z` with single rounding error, as defined by the platform C library. The result is the same as if the operation was carried with infinite precision and rounded to a floating-point number.'''

Attached is a simple module for Python 3 demonstrating the fused multiply-add operation. On my machine, example.py prints:

40037.524591982365   horner_double
40037.48821639768    horner_fma
40037.49486325783    horner_compensated
40037.49486325783    horner_decimal

[1] https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation
[2] S. Graillat, P. Langlois, N. Louvet. Accurate dot products with FMA. 2006
[3] S. Graillat, Accurate Floating Point Product and Exponentiation. 2007.
[4] S. Graillat, P. Langlois, N. Louvet. Improving the compensated Horner scheme with a Fused Multiply and Add. 2006
[5] J.-M. Muller, N. Brisebarre, F. de Dinechin, C.-P. Jeannerod, V. Lefèvre, G. Melquiond, N. Revol, D. Stehlé, S. Torres. Handbook of Floating-Point Arithmetic. 2010. Chapter 5
[6] ISO/IEC 9899:TC3, "7.12.13.1 The fma functions", Committee Draft - Septermber 7, 2007
[7] https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
[8] https://www.python.org/dev/peps/pep-0007/

Thread on python-ideas:
https://mail.python.org/pipermail/python-ideas/2017-January/044300.html

What's the point of adding this in the math module rather than a more specialized library like Numpy?

I would say because it has wide applicability, especially considering the amount of code it adds. It is similar in spirit to copysign or fsum which are already there and C99 includes it anyway. For simpler things like dot product or polynomial evaluation importing Numpy might be too much of an dependency.

I don't know. If I want to compute a dot product, the first thing I'll do is import Numpy and then use the @ operator on Numpy arrays.

The performance argument unlikely is applicable in this case. I suppose that an overhead of function call in Python make two operators faster than one function call.

Alternatives to fma() for exact computations are integer arithmetic (if all values can be scaled to integers), fractions and decimal numbers.

But since fma() is a part of C (C99), C++ (C++11) and POSIX (POSIX.1-2001) standards for long time, I don't have objections against including it in the math module.

The performance argument unlikely is applicable in this case.

Agreed. This is mainly about accuracy, not speed: the FMA operation is a fundamental building block for floating-point arithmetic, is useful in some numerical algorithms, and essential in others (especially when doing things like double-double arithmetic). It would be valuable to have when prototyping numerical algorithms in pure Python.

Given that it's supported in C99 and on current Windows, I'm +1 on including it in the math module.

Note that implementing this it not quite as straightforward as simply wrapping the libm version, since we'll also want the correct exceptional behaviour, for consistency with the rest of the math module: i.e., we should be raising ValueError where the fma operation would signal the invalid FPE, and OverflowError where the fma operation would signal the overflow FPE.

An implementation note: IEEE 754-2008 leaves it to the implementation to decide whether FMA operations like:

0 * inf + nan

and

inf * 0 + nan

(where nan represents a quiet NaN and the inf and 0 can have arbitrary signs) signal the invalid operation FPE or not. (Ref: 7.2(c) in IEEE 754-2008.)

I'd suggest that in this case we follow what Intel does in its x64 chips with FMA3 support.(according to ). If I'm reading the table in section 2.3 of the Intel Advanced Vector Extensions Programming Reference correctly, Intel does *not* signal the invalid operation FPE in this case. That is, we're following the usual rule of quiet NaN in => quiet NaN out with no exception. This does unfortunately conflict with the IBM decimal specification and Python's decimal module, where these operations *do* set the "invalid" flag (see the spec, and test fmax0809 in the decimal test set).

Isn't the behaviour of quiet NaNs kindof implementation-dependent already?

Isn't the behaviour of quiet NaNs kindof implementation-dependent already?

Not as far as IEEE 754-2008 is concerned, and not as far as Python's math module is concerned, either: handling of special cases is, as far as I know, both consistent across platforms and compliant with IEEE 754. That's not true across Python as a whole, but it should be true for the math module.

If you find an exception to the above statement, please do open a bug report!

Here's a pure Python reference implementation, with tests.

And here's a patch.

FYI :
dotnet/runtime#98704

It's worth noting that if you want a performant and open source implementation, you have:

The official AMD's implementation: https://github.com/amd/aocl-libm-ose/blob/master/src/isa/avx/fmaf.c
- In plain C, portable, no intrinsics
- BSD 3-Clause
The Berkeley SoftFloat implementation: https://github.com/ucb-bar/berkeley-softfloat-3/blob/master/source/s_mulAddF32.c
- In plain C, portable, no intrinsics
- Not as well optimized

I believe there is also one available from netlib, but I don't have that link handy. Arm typically has algorithms available on its repository, but doesn't provide one for fma because all Arm64 based hardware provide an fma instruction natively: https://github.com/ARM-software/optimized-routines

@tannergooding

@tannergooding double(64bits float) version is desired, maybe :
https://github.com/amd/aocl-libm-ose/blob/master/src/isa/avx/fma.c
https://github.com/ucb-bar/berkeley-softfloat-3/blob/master/source/s_mulAddF64.c

Yes, if you want the 64-bit (double) versions they exist in neighboring files. Notably doing something like (float)fma((double)x, (double)y, (double)z) may not be sufficient in all cases if you want a 32-bit (float) result. There is a subset of operations for which the rounding twice is known to be a non-issue, but I don't recall off the top of my head if fma is included in that, so it may be desirable to have both (and there may even be a performance difference).

pythongh-73468: Add math.fma() function

My previous attempt in 2020 failed. I try again in 2024: PR #116667

pythongh-73468: Add math.fma() function

gh-73468: Add math.fma() function (#116667)

Function added by commit 8e3c953. I only took 7 years only Linux, macOS and Windows get a good support of fma() :-)

pythongh-73468: Add math.fma() function (python#116667)

Navigation Menu

Repository navigation

Fused multiply-add: proposal to add math.fma() #73468

Linked PRs

27 remaining items

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Fused multiply-add: proposal to add math.fma() #73468

Description

Linked PRs

Activity

jurajsukop commented on Jan 16, 2017

vstinner commented on Jan 16, 2017

pitrou commented on Jan 16, 2017

jurajsukop commented on Jan 16, 2017

pitrou commented on Jan 16, 2017

serhiy-storchaka commented on Jan 16, 2017

mdickinson commented on Jan 16, 2017

mdickinson commented on Jan 16, 2017

pitrou commented on Jan 16, 2017

mdickinson commented on Jan 16, 2017

mdickinson commented on Jan 16, 2017

mdickinson commented on Jan 17, 2017

27 remaining items

clairvoyante commented on Feb 21, 2024

tannergooding commented on Feb 22, 2024

clairvoyante commented on Feb 23, 2024

tannergooding commented on Feb 23, 2024

vstinner commented on Mar 13, 2024

vstinner commented on Mar 17, 2024

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions