[Python] How to do fft based convolution with long doubles/float128s

On my linux system I have:

np.finfo(np.float128)
info(resolution=1e-18, min=-1.189731495357231765e+4932, max=1.189731495357231765e+4932, dtype=float128)


So that is an 80-bit long double. I want to perform a convolution between two reasonably long arrays of np.float128s. scipy.signal.convolve works with method='direct' but gives the wrong answer for method='fft'. Here is a toy example:

a = np.array(['1.e+401', '1.e+000', '1.e+401', '1.e+000'], dtype=np.float128)

convolve(a, a, mode='full', method='direct')
array([1.e+802, 2.e+401, 2.e+802, 4.e+401, 1.e+802, 2.e+401, 1.e+000],
dtype=float128) # correct
convolve(a, a, mode='full', method='fft')
array([1.e+802, 0.e+000, 2.e+802, 0.e+000, 1.e+802, 0.e+000, 0.e+000],
dtype=float128) # wrong


I tried implementing the convolution from scratch using pyfftw but it still gave the wrong answer.

How can this be done?

Continue reading...