Skip to content

FFT with kissfft backend does not work correctly with extended precision floating numbers

Hi!

The implementation of kissfft is uses "double" in one line, but is otherwise scalar type agnostic. https://gitlab.com/libeigen/eigen/-/blob/c29c800126982c561e8d0b9255dc65474cd98de3/unsupported/Eigen/src/FFT/ei_kissfft_impl.h#L35

And of course this matters where double does not have enough precision, say if I try to use arbitrary precision arithmetic as suggested for example here: https://scicomp.stackexchange.com/questions/23867/c-libraries-for-fast-fourier-transform-in-high-precision I didn't really check, but perhaps this gives wrong results already for computations with long double.

There's a trivial fix:

diff --git a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
index 430953aee..f1beeb227 100644
--- a/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
+++ b/unsupported/Eigen/src/FFT/ei_kissfft_impl.h
@@ -31,7 +31,7 @@ struct kiss_cpx_fft
     using numext::cos;
     m_inverse = inverse;
     m_twiddles.resize(nfft);
-    double phinc =  0.25 * double(EIGEN_PI) / nfft;
+    Scalar phinc = numext::atan(Scalar(1)) / nfft;
     Scalar flip = inverse ? Scalar(1) : Scalar(-1);
     m_twiddles[0] = Complex(Scalar(1), Scalar(0));
     if ((nfft&1)==0)

I'm not sure this is the right way to do things, because one has to compute PI on each call to FFT. Although, this is just one call to atan, and kissfft is not competing for speed from what I understand, so maybe it is OK as a fix?