FFTW und doppelter Präzision

Question

ich die FFTW 3.3.5 Bibliotheken mit (http://www.fftw.org/doc/Precision.html) zusammengestellt:

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

ich kompilieren Sie den Code unten mit 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])); 
} 

Um meine das Ergebnis der endgültigen Verständnis printf muss die gleichen bis zu 17-18 Dezimalziffern für alle Läufe sein, aber was ich bekomme ist in der 14. Dezimalstelle deutlich. So etwas könnte bedeuten, dass ein langer doppelter zu einem doppelten Typ degradiert wurde. Die Ausgabe der Code ändert sich von einem Lauf zum anderen:

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 

Irgendwelche Ideen auf, wo ich verloren die lange Doppel Genauigkeit?

Answer 1

Jedes Element des Arrays, das aus einer komplexen FFT resultiert, ist a weighted sum of the 4096 items of the input array, das auf verschiedene Arten berechnet werden kann. Eine Untersuchung der Art, wie sich der Fehler ausbreitet, kann wie beschrieben here durchgeführt werden.

Da die Gewichte der Summe immer der Norm 1 entsprechen, ist die Varianz eines Elements des Ausgabearrays die Summe der Varianz der Elemente des Eingabearrays (Elemente des Eingabearrays werden als unabhängige Zufallsvariablen betrachtet). .

Die Präzision 1e-18 ist das Verhältnis der Standardabweichung zu der Norm des Wertes. Daher werden, wenn lange Doppel verwendet:

Schließlich ist die Genauigkeit des Ergebnisses schreibt:

Eine direkte Folge dieser Gleichung das Phänomen der katastrophalen Kündigung ist: der Standard Abweichung der Ausgabe ist endlich, aber der Wert |y| ist sehr klein, da ähnliche Werte subtrahiert wurden ... Folglich ist keine der Ziffern signifikant! Siehe hier (https://en.wikipedia.org/wiki/Loss_of_significance) zum Beispiel. Dies erklärt, warum 7.463682637972633967E-24 in Ihrem Fall leicht 1.965167197605865111E-23 werden kann. Tatsächlich kann, wenn eine teilweise Aufhebung auftritt, die Genauigkeit leicht von 1e-18 auf Werte wie 1e-14 abnehmen. Kleinere Werte von | y | weisen weniger signifikante Ziffern auf.

Die Tatsache, dass sich die Ausgabe von einem Lauf zum anderen unterscheidet, kann auf die Verwendung des Flags FFTW_MEASURE zurückzuführen sein. Wenn dieses Flag verwendet wird, werden verschiedene Algorithmen ausprobiert und FFTW wählt das schnellste aus. Da die 1D-FFT der Länge 4096 nun sehr schnell ist, ist das Timing möglicherweise nicht sehr konsistent und verschiedene Algorithmen vergleichbarer Leistungen können ausgewählt werden. Die verschiedenen Algorithmen führen zu unterschiedlichen Berechnungen und unterschiedlichen Ergebnissen für nicht signifikante Ziffern. Bestehen diese Varianten noch, wenn stattdessen das Flag verwendet wird? Wenn dieses Flag verwendet wird, it seems that FFTW uses a simple heuristic to choose the algorithm. Es ist wahrscheinlich, dass diese Heuristik deterministisch ist ...