#include typedef double _Complex c; typedef long double _Complex cx; void fft(const int n, c a[n]){ {int j=0; for(int i=0; i #define LN 10 #define N (1< void pc(c a){printf("%14.7e + %14.7ei\n", creal(a), cimag(a));} int main(){c w[N]; // Initialize array w, thereafter a constant. {c W[LN]; {W[LN-1] = -1; for(int j=LN-1; j; --j) W[j-1] = csqrtl(W[j]);} {int j=N; while(j--) {int k=LN; c a = 1; while(k--) if(j & (1 << k)) a *= W[k]; w[j] = a;}}} // Demonstrate that w[1] is the correct root of unity // and then that the other values of w are correct. {c r = w[1]; printf("r = "); pc(r); {int j = N; while(j--) { c e = w[j]*r - w[(j+1)&(N-1)]; if(fabs(__real__ e) > ep || fabs(__imag__ e) > ep) { pc(w[j]); pc(w[(j+1)&(N-1)]); printf("Roots: %d\n", j);}}} {int j = LN; while(j--) r *= r; printf("r^%d-1 = ", N); pc(r-1);} } // Evaluate transform for each basis vector and // demonstrate that the answer is correct. {int k=N; while(k--) {c A[N]; int j=N; while(j--) A[j]=0; A[k]=1; fft(N, A); {int n=N; while(n--) { c p = w[n*k &(N-1)]; if(fabs(creal(A[n]) - creal(p)) > ep || fabs(cimag(A[n]) - cimag(p)) > ep) {pc(A[n]); pc(p); printf("er: %6d %6d\n", k, n);}}}}} return 0;} /* r = 9.9998118e-01 + 6.1358846e-03i r^1024-1 = -2.5313085e-14 + -6.6613381e-16i */