FFTW和长双精度

我（http://www.fftw.org/doc/Precision.html）编制的FFTW 3.3.5库：FFTW和长双精度

./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位小数不同。像这样的东西可能表示长双被降为双类型。从一个运行代码切换到另一个输出：

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

上，我失去了长双精度任何想法？

来源

2016-12-10 RunsWithScissors

由复数FFT产生的数组中的每一项都是a weighted sum of the 4096 items of the input array，可以用不同的方式计算。对错误传播方式的研究可以按照here所述进行。

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

$\Delta y^2=\sum\Delta x^2$

精度1E-18是标准偏差的值设定的范数的比率。因此，如果长期双重使用：

$\Delta y^2=\sum10^{-36}\left|x\right|^2$

最后结果的精度写道：

$\frac{\Delta y}{|y|}=10^{-18}\frac1{|y|}\sqrt{\sum\left|x\right|^2}$

这个方程的一个直接后果是灾难性的取消的现象：标准输出的偏差是有限的，但值|y|是非常小的，因为类似的值被减少了...因此，没有一个数字是显着的！例如，请参阅（https://en.wikipedia.org/wiki/Loss_of_significance）。这就解释了为什么7.463682637972633967E-24可以很容易地变成你的情况下的1.965167197605865111E-23。事实上，如果发生部分取消，精度可以从1e-18轻松降至1e-14等值。较小的| y |值具有较少有效位数。

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

来源

2016-12-26 00:14:03 francis

FFTW和长双精度

回答

相关问题