[SciPy-Dev] Accuracy of single-precision FFT

[Sturla Molden]

Den 25.06.2010 16:05, skrev Pauli Virtanen:
>
> Supporting only MKL is not an option, since even if we managed to secure
> ourselves a redistribution license,
>    


I was not suggesting that :) Or at least I did not intend to.

We have to have other alternatives. For example to run NumPy on ARM. Or 
if someone wants to py2exe a NumPy based program, and does not own an 
MKL license themselves.

But for example: If FFTPACK is faulty on single precision FFTs, we could 
restrict its use to double precision. MKL could still be used for single 
precision.

I also wonder if this in-accuracy is due to FFTPACK's algorithms or the 
Fortran compiler? NumPy uses  fftpack_lite which is a C library. It is 
easy to compile fftpack_lite for single precision as well (basically 
just tydef'ing Treal to float). Is that inaccurate too? If not, we have 
a Fortran compiler problem.

With C99 or C++, we can even compile fftpack_lite for complex numbers, 
it's just a litte bit more work than:

typedef double _Complex Treal; // C99
typedef std::complex<double> Treal;  // C++

What's left then is the DCT. I'm sure we can get around that as well. We 
can e.g. compute DCT via an FFT.

This means Scipy does not need FFTPACK at all anymore, since all FFTPACK 
does is in fftpack_lite and the C99 or C++ compiler. All in all, 
FFTPACK  is close to obsolete.

Note that NumPy's fftpack_lite has a better wrapper as well. The FFTPACK 
wrapper leaks memory (for plans) and prevents the GIL from being 
released (the plan cache is not thread safe). With fftpack_lite however, 
we can safely release the GIL (I've submitted code for that before), as 
the plans are cached in Python and not C. With fftpack_lite we can also 
empty the plan cache by rebinding the dict to an empty one.

Note that with fftpack_lite it is also trivial to compile the FFT for 
long double. We have those dtypes in NumPy too.


Sturla