如何快速，准确地使用长双打/float128进行基于FFT的卷积

Question

在我的Linux系统中，我有：

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

这是一个80位长的双倍。我想在两个合理的长阵列之间进行卷积。

np.float128

与

scipy.signal.convolve

一起工作，但给出了错误的答案。这是一个玩具示例：

method='direct'

我尝试使用

method='fft'

从头开始实施卷积，但它仍然给出了错误的答案。 我的最终目标是能够使用FFT快速准确地进行类似的卷积：

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

怎么做？

Scipy似乎实际上并不支持

pyfftw

数据类型，因为在适当支持此类型的Linux计算机上使用该类型时，我会遇到以下非常清晰的错误。 eDeprecationWarning：dtype = float128不受相关性的支持，并且会引起Scipy 1.17.0的错误。支持的dtypes是：布尔，整数，

a = np.array([1e401, 1e000, 1e401, 1e000] * 10000, dtype=np.float128)
convolve(a, a)

，

Answer 1

，

np.float16

，

np.float32

，

np.float64

我的Scipy版本是1.15.1，它是撰写时PYPI

上最latat的版本。

技术上，Pyfftw应该支持它，因为该文档明确指出： 在内部，支持FFT的三个精度。这些对应于单个精度浮点，双精度浮点和长双精度浮点，这对应于Numpy的Float32，Float64和

longdouble

dtypes（和相应的复杂类型）。精度由相关方案决定，该方案由输入数组的DTYPE指定。 事情是，这不是一个精确的问题，而是数值稳定性问题

！确实，让我们检查输出数量较小的双精度数字：

np.complex64 详细输出是：

np.complex128

这表明当指数很大时，FFT方法不再正确。实际上，尽管具有双精度数据类型的输入值，数值误差已经很大，但诸如

a = np.array([1e40, 1, 1e40, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Wrong
print()

a = np.array([1e5, 1, 1e5, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Correct
print()

a = np.array([1e15, 1, 1e15, 1])
print(convolve(a, a, mode='full', method='direct'))   # Correct
print(convolve(a, a, mode='full', method='fft'))      # Wrong, but... Oh!

的输入值支持高达1E+308的值。实际上，即使输入值的

[1.e+80 2.e+40 2.e+80 4.e+40 1.e+80 2.e+40 1.e+00]
[1.e+80 0.e+00 2.e+80 0.e+00 1.e+80 0.e+00 0.e+00]

[1.e+10 2.e+05 2.e+10 4.e+05 1.e+10 2.e+05 1.e+00]
[1.e+10 2.e+05 2.e+10 4.e+05 1.e+10 2.e+05 1.e+00]

[1.e+30 2.e+15 2.e+30 4.e+15 1.e+30 2.e+15 1.e+00]
[1.00000000e+30 1.99520395e+15 2.00000000e+30 3.93577438e+15 1.00000000e+30 1.94544573e+15 0.00000000e+00]

输入值与双重精度支持相比，这是错误的！当指数太大（> 〜7.97）时，某些值会夹紧到0，结果显然不再可靠。

基于FFT的算法在数值上是不准确的。这很可能是由于催化性的取消

由于在算法中多次被小值减去巨大的值。提高精度有所帮助，但是在80位x86特异性

1e15

数据类型中可用的少数数字远非在这里足够（例如64位的数据类型远非足以使输入值是1E15）。

实际上，如果我们使用

1e8

数据类型，我们实际上只能看到一个略有改进：

long double

我认为当问题大致出现时

np.float128

。 我不知道是否可以解决此数值稳定性问题，还是这是使用FFT进行卷积的基本问题。话虽这么说，请注意，您可以通过调用支持

Quad-precision编号的本机FFTW实现来进一步提高精度（在带有GCC的Linux上，也可能也可能是Clang）。这应该有助于进一步推动极限：数字增加两倍。这一额外的精度以高昂的价格出现：较慢的计算速度较慢，因为四元素数字不是本地数字。实际上，在X86-64 CPU上（甚至在其他体系结构上不支持），80位浮点（FP）数量已经很慢，因为计算是在X87 FP单元上执行的，该单元自1〜20年以来就相当过时了. 避免此数值问题的密钥可能是避免输入数字的巨大变化。这是afaik称为先进（很常见）。一种常见的做法是计算日志空间中的值，但是在这种情况下，这不再是卷积...因此，我不知道在这种情况下是否可以进行预处理。