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和長雙精度

回答

相關問題