math.gamma result slightly different on aarch64-apple-darwin

[youknowone]

Bug report

Bug description:

The difference is watched from distributions.
Because the result is changed by fp-contract. It may not be strictly a bug in source code, but can be a bug of build & distrubition.

On x86_64 linux:

Python 3.13.3 (main, Apr 21 2025, 11:31:07) [GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import math; math.gamma(-3.8510064710745118)
0.36222232384328096

On aarch64 macOS:

Python 3.13.2 (main, Feb  4 2025, 14:51:09) [Clang 16.0.0 (clang-1600.0.26.6)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import math; math.gamma(-3.8510064710745118)
0.36222232384328107

tgamma.c file is separated from mathmodule.c to test easy.
new tgamma.c

The result below is the result of tgamma.c, but they are the same in CPython distribution.
Test results:

tgamma.c in x86_64 linux:

$ clang -lm tgamma.c -ffp-contract=off && ./a.out The result of tgamma(-3.8510064710745118) is: 0.362222
Bit representation of result: 4600196837208649421
$ clang -lm tgamma.c -ffp-contract=on && ./a.out 
The result of tgamma(-3.8510064710745118) is: 0.362222
Bit representation of result: 4600196837208649421

tgamma.c in aarch64 macos:

$  clang -lm tgamma.c -ffp-contract=off && ./a.out
The result of tgamma(-3.8510064710745118) is: 0.362222
Bit representation of result: 4600196837208649421
$  clang -lm tgamma.c -ffp-contract=on && ./a.out
The result of tgamma(-3.8510064710745118) is: 0.362222
Bit representation of result: 4600196837208649423

The original test I found this mismatch:
On aarch64-apple-darwin,

x is input.
Left is CPython, right is pymath

thread 'gamma::tests::test_tgamma' panicked at src/gamma.rs:303:17:
assertion `left == right` failed: x = -3.8510064710745118, py_gamma = 0.36222232384328107, rs_gamma = 0.36222232384328096
  left: 4600196837208649423
 right: 4600196837208649421

thread 'gamma::tests::test_tgamma' panicked at src/gamma.rs:303:17:
assertion `left == right` failed: x = -3.8510064710745118, py_gamma = 0.36222232384328107, rs_gamma = 0.36222232384328096
  left: 4600196837208649423
 right: 4600196837208649421

thread 'gamma::tests::test_tgamma' panicked at src/gamma.rs:286:5:
Test failed: assertion `left == right` failed: x = -3.8510064710745118, py_gamma = 0.36222232384328107, rs_gamma = 0.36222232384328096
  left: 4600196837[208](https://github.com/RustPython/pymath/actions/runs/14565988158/job/40855284918?pr=2#step:6:209)649423
 right: 4600196837208649421.
minimal failing input: x = -3.8510064710745118
	successes: 0
	local rejects: 0
	global rejects: 0

On x86_64-unknown-linux-gnu,

x = -3.8510064710745118
py_gamma = 0.36222232384328096
rs_gamma = 0.36222232384328096
py_gamma_repr = 4600196837208649421
rs_gamma_repr = 4600196837208649421

Turning off fp-contract will fix aarch64 macos build to match x86 linux/windows build.

CPython versions tested on:

3.13

Operating systems tested on:

macOS

[corona10]

cc @rhettinger @ned-deily

[tim-one]

It's a 2 ULP difference:

>>> 0.36222232384328096 .hex()
'0x1.72ea68ab26ecdp-2'
>>> 0.36222232384328107 .hex()
'0x1.72ea68ab26ecfp-2'

I don't consider this to be "a bug". CPython makes no promises about the accuracy of "advanced" math functions,, and is typically at the mercy of the platform's C compiler and libm. They can't be expected to deliver identical results across platforms, compilers, or even across releases (of the compiler or of Python) on a single platform.

In this case, looks like the MacOS C compiler is exploiting the HW's fused multiply-add (FMA) support. There's nothing in our source code that asks for that, or that forbids it. It's up to the compiler, and that's fine. FMA is generally a good thing, running faster and with smaller error. Of course it can have tiny effects on the low-order bits of the final result (or, for algorithms that rely on it, which Python doesn't use, huge effects).

[skirpichev]

FMA is generally a good thing, running faster and with smaller error.

BTW, in the OP case - aarch64 result seems slightly better (rather 1ULP off).

The difference is watched from distributions.

It depends on how you do this. If you compare computed value with some float by == - that's rather bug in your tests.

[youknowone]

I don't consider this to be "a bug".

I agree this is not necessarily a bug. Checking whether this is a policy of CPython will be good enough. But it is not because it is 2ULP difference, but because no promise about advanced math funcitions, right?

Finding bigger difference is possible.

$ clang -lm tgamma.c && ./a.out
The result of tgamma(-3.8510064710745118) is: 3.809507
Bit representation of result: 4615760666384231000

$ clang -lm tgamma.c && ./a.out
The result of tgamma(3.6215752811868267) is: 3.809507
Bit representation of result: 4615760666384231005

If you compare computed value with some float by ==

I am not comparing the float values but comparing the bit representation of float values. Probably I didn't write the issue detailed enough. I am not insisting the result of math.gamma is incorrect, but only saying the bit representation of gamma function is different only in aarch64-apple-darwin unlike x86 linux or windows on major platforms I have.

For more context, I was going to write the port of those functions in Rust. I know float operation result can be varied by based platform, libm or a few other reasons.
The rust implementation, using same libc functions and carefully ported code produced the same bit representation in all of the platforms. The bit representation was also same on CPython on x86-64 linux and windows. Surprisingly, only CPython on aarch64-apple-darwin had different bit representation to the other 3 rust + 2 cpython results. I still think this is not strictly a bug. So I'd like to check if this is intended or not intended. I thought there could be a chance a build flag was missed for aarch64.

[skirpichev]

Finding bigger difference is possible

Sorry, it seems you use different arguments in your latest examples. -3.8510064710745118 vs 3.6215752811868267. It's hard to say in which example the difference is bigger.

but comparing the bit representation of float values.

It's same story. Rather you should use something like math.isclose().

[youknowone]

Sorry, the log is hardcoded one. (The source code is linked above)
3.6215752811868267 is the correct input in the comment. You can see bit representations of them to check difference.

@skirpichev I am not getting why you are emphasizing math.isclose() in same story.
I know comparing bits is not the best practice for float operations. I described what was the purpose of the issue and I am not insisting the difference of small values are real bugs. If it was not clear enough, please tell me which part was misleading. I will try to fix it.

[corona10]

Or leave note at docs.python.org that the it could show different bit results according to the compiler options?
Because we can not control all distros that are using custom compiler options.

[youknowone]

@corona10 Different bit representations for floating-point operations across platforms are expected behavior in computer science. So I feel like that's a redundant commentary for general computer science knowledge.

If CPython wants to try best to keep consistent cross-platform floating-point behavior, modifying the build configuration with appropriate compiler flags would be the solution.

However, if the current default compiler options align with the project's policy of minimal compiler customization, then I believe we can simply close this issue as working as designed.

[skirpichev]

3.6215752811868267 is the correct input in the comment. You can see bit representations of them to check difference.

Thanks, I guess in the second case you have result on aarch64 (on x86_64 linux box I got your first result). In this case, this is worse by 1ULP (3ULP vs 2ULP) than x86_64 result:

In [17]: struct.unpack('<d', struct.pack('<q', 4615760666384231000))  # 1st example
Out[17]: (3.8095071915716225,)

In [18]: struct.unpack('<d', struct.pack('<q', 4615760666384231005))  # 2nd
Out[18]: (3.8095071915716248,)

In [19]: f"{gamma(3.6215752811868267):.17f}"  # mpmath, computed with mp.prec=113
Out[19]: '3.80950719157162333'

In [20]: 3.80950719157162333-3.8095071915716225
Out[20]: 8.30000000000000000061569678209333068e-16

In [21]: 3.80950719157162333-3.8095071915716248
Out[21]: -1.47000000000000000001158839675115145e-15

In [22]: math.ulp(3.80950719157162333)
Out[22]: 4.440892098500626e-16

In [23]: _20/_22
Out[23]: 1.86899384535875584013864258992043688

In [24]: _21/_22
Out[24]: -3.31014572611731456002609474964515357

In [25]: math.gamma(3.6215752811868267)  # my x86_64 box (c.f. In[17])
Out[25]: 3.8095071915716225

BTW, could you test libm's tgamma() on your systems?

why you are emphasizing math.isclose() in same story.

It's a good idea to use something like this, unless you are sure that your computations are identical (i.e differ only by transformations, that are true for all floats, e.g. you can't use associativity to reorder operands). Optimization, using HW's FMA - is an example of transformation (permitted with -ffp-contract=on), which can break this.

modifying the build configuration with appropriate compiler flags would be the solution.

I doubt that forbidding floating-point contraction is a good idea. Besides a performance boost you will, probably, get more accurate results on average.

[tim-one]

[@youknowone wrote]

Checking whether this is a policy of CPython will be good enough. But it is not because it is 2ULP difference, but because no promise about advanced math funcitions, right?

Right! Although, to be fair, CPython makes no guarantees about the results of simple float arithmetic (+ - * /) either. Those results are also inherited from the platform C. It just so happens that almost all platforms today have HW that implements the IEEE-754 standard, which does define the precise results of simple arithmetic.

gamma is a little different in that CPython implements it itself. Our code doesn't rely on FMA, but it's up to the compiler whether to use it in lines like:

            num = num * x + lanczos_num_coeffs[i];
            den = den * x + lanczos_den_coeffs[i];

Results can differ.

Why does CPython implement gamma itself? Because some platforms' implementations at the time were too sloppy to bear. From the comments:

/* Implementation of the real gamma function.  Kept here to work around
   issues (see e.g. gh-70309) with quality of libm's tgamma/lgamma implementations
   on various platforms (Windows, MacOS).  In extensive but non-exhaustive
   random tests, this function proved accurate to within <= 10 ulps across the
   entire float domain.  Note that accuracy may depend on the quality of the
   system math functions, the pow function in particular.  Special cases
   follow C99 annex F.  The parameters and method are tailored to platforms
   whose double format is the IEEE 754 binary64 format. ...
*/

Note that errors as large as 10 ULP were seen at the time it was checked in. It's possible that the platforms in question have improved their implementations since then. The kind of testing done in gh-70309 could be done again to see whether we still "need" to supply our own.

[youknowone]

Thank you all. I am understanding the situation better.

        num = num * x + lanczos_num_coeffs[i];
        den = den * x + lanczos_den_coeffs[i];

I confirmed this is the exact line why this happned. When I changed this operation to fma, tgamma always returns the "different" value of aarch64 macos. lgamma seems to have more place than that.

I also confirmed this is only happening on aarch64 macos, not on x86_64 macos. Probably compiler default is set to be more aggressive for aarch64.

Rather than changing compiler ops, how do you think making some operations in gamma to be explicitly using fma or hostile to be automatically fma applied by which one is more accurate to real function? Explicitly using fma will make FMA-supporting architectures to return always the same value and another group works as before. Making hostile to fma will make everything agree on the other way.
Of course nothing need to be done if it is not worth to do.