rfft output difference python 2.7.15 / python 3.7.2

[julienprieur]

Same code using rfft produces slighly different output on python 2 / python 3.
I noticed it because it broke some automatic testing of my internal libraries.

It looks like the latest digit of the real part is being truncated in python 3.

Reproducing code example:

import numpy as np
print(np.fft.rfft(np.sin(np.arange(1000) * 0.1))[25])

output:
(-0.686231256675454-5.4138146836022445j) in python 3.7.2
(-0.6862312566754544-5.4138146836022445j) in python 2.7.15

Numpy/Python version information:

Numpy version 1.16.0
Python 2.7.15 x64
Python 3.7.2 x64
Windows 10

[chrisjbillington]

If you're comparing the numbers based on how they print, it could just be a change in how floats are printed between Python 2 and Python 3, see here (emphasis mine):

Interestingly, there are many different decimal numbers that share the same nearest approximate binary fraction. For example, the numbers 0.1 and 0.10000000000000001 and 0.1000000000000000055511151231257827021181583404541015625 are all approximated by 3602879701896397 / 2 ** 55. Since all of these decimal values share the same approximation, any one of them could be displayed while still preserving the invariant eval(repr(x)) == x.

Historically, the Python prompt and built-in repr() function would choose the one with 17 significant digits, 0.10000000000000001. Starting with Python 3.1, Python (on most systems) is now able to choose the shortest of these and simply display 0.1.

Note that this is in the very nature of binary floating-point: this is not a bug in Python, and it is not a bug in your code either. You’ll see the same kind of thing in all languages that support your hardware’s floating-point arithmetic (although some languages may not display the difference by default, or in all output modes).

If you need the numbers to be identical across Python versions, you should compare their binary representations to see if they are equal, not their string representations. (you might be doing this already and the printing is just for the bug report, but for what it's worth)

Then again, I see a different result again on my setup, and I see this both on Python 2 and Python 3:

>>> print(np.fft.rfft(np.sin(np.arange(1000) * 0.1))[25])
(-0.6862312566754543-5.413814683602245j)

But I am using the intel math kernel library, perhaps it has a different implementation of sin, or a different fft implementation. You would need to double check that both your Python 2 and 3 installations are using the exact same underlying math libraries. Whilst there are standards defining what basic arithmetic should do for floats, I don't think the same thing exists for ffts or the sine function, such that different implementations may have slightly different floating point rounding errors.

[julienprieur]

I'm using the stock python and numpy verison (not using MKL).

I did some further tests and it turned out that the difference likely comes from different underlying math libraries. I tested the sin function from the python math package and I get differences on latest digit for about 2% of inputs.
I guess that the difference is due to the fact that both python versions are linked against different c runtime (python 2.7.15 is compiled with MSC v1500 vs v1916 for python 3.7.2).

I think this issue can be closed then as not linked to numpy issue.

[vkhodygo]

@chrisjbillington AFAIK, Intel indeed uses its own implementation of fft and that is clearly understandable. However, I'm not sure about basic math functions. I get exactly the same answer (icc build as well).

@julienprieur there is a very nice answer that explains some differences. I also have checked stock numpy (using a VM) and got the same results.