FFTW 和 long 双精度答案

【问题标题】：FFTW and long double precisionFFTW 和 long 双精度
【发布时间】：2016-12-10 16:58:56
【问题描述】：

我使用 (http://www.fftw.org/doc/Precision.html) 编译了 FFTW 3.3.5 库：

./configure --enable-long-double
make
make install

我用gcc -std=gnu99 main.c -o sample.x -lfftw3l -lm编译下面的代码

#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include <string.h>

#define PI acosl(-1.0L)
#define FMODE FFTW_MEASURE

int main() {
  fftwl_complex  *A = fftwl_malloc(4096*sizeof(fftwl_complex));
  fftwl_plan     ft = fftwl_plan_dft_1d(4096, A, A, FFTW_BACKWARD, FMODE);
  long double    q, u, overN = ((long double) 1.L/4096);

  for (long int j = 0; j < 4096; j++) {
    q = 2.L*PI*(j*overN - 0.5L);
    u = 2.L*atan2l(0.5L*sinl(0.5L*q),cosl(0.5L*q));
    A[j] = -1.IL*cpowl(0.01L*(1.L/ctanl(0.5L*(u-0.1IL)) - 1.IL),2);
  }
  printf("%26.18LE\t%26.18LE\n", creall(A[1]), cimagl(A[1]));
  fftwl_execute(ft);
  for (int j = 0; j < 2048; j++) {
    A[j] = -1.0IL*((fftwl_complex) j*A[j])*overN;
  }
  printf("%26.18LE\t%26.18LE\n", creall(A[1]), cimagl(A[1]));
  memset(A+2048, 0, 2048*sizeof(fftwl_complex));
  fftwl_execute(ft);
  printf("%26.18LE\t%26.18LE\n", creall(A[1]), cimagl(A[1]));
}

据我了解，最终 printf 的结果必须是但是，所有运行都使用相同的最多 17-18 位十进制数字我得到的在小数点后 14 位是不同的。这样的事情可能表明 long double 被降级为 double 类型。代码的输出从一次运行到另一次运行的变化：

  2.907416794556517046E-07    9.025765251354815022E-05
 -5.697284273172913999E-04    7.463682637972633967E-24
  1.895341327532694420E-04    3.343168537700265992E-07

  2.907416794556517046E-07    9.025765251354815022E-05
 -5.697284273172913999E-04    1.965167197605865111E-23
  1.895341327532697672E-04    3.343168537692799396E-07

关于我失去long double准确性的任何想法？

【问题讨论】：

标签： c

【解决方案1】：

从复数 FFT 得到的数组中的每一项都是a weighted sum of the 4096 items of the input array，可以用不同的方式计算。可以按照here 的描述来研究错误传播的方式。

由于总和的权重始终为范数 1，因此输出数组的项的方差是输入数组项的方差之和（输入数组的项被视为独立随机变量）。

精度 1e-18 是标准偏差与值的范数之比。因此，如果使用 long double：

最后，精度写成：

这个等式的直接后果是灾难性抵消现象：输出的标准偏差是有限的，但值|y| 非常小，因为减去了相似的值......结果，没有一个数字很重要！例如，请参见此处（https://en.wikipedia.org/wiki/Loss_of_significance）。这解释了为什么 7.463682637972633967E-24 在您的情况下很容易变成 1.965167197605865111E-23。实际上，如果发生部分抵消，精度很容易从 1e-18 降低到 1e-14 之类的值。 |y| 的较小值具有较低的有效数字。

每次运行的输出不同的事实可能是由于使用了标志FFTW_MEASURE。如果使用此标志，将尝试不同的算法，FFTW 将选择最快的算法。由于长度为 4096 的一维 FFT 现在非常快，因此时序可能不是很一致，并且可以选择性能相当的不同算法。对于非有效数字，不同的算法将导致不同的计算和不同的结果。如果改用标志FFTW_ESTIMATE，这些变化是否仍然存在？如果使用此标志，it seems that FFTW uses a simple heuristic to choose the algorithm。这种启发式方法很可能是确定性的……

【讨论】：