7 static void make_sintbl(int n, float sintbl[])
10 double c, s, dc, ds, t;
12 n2 = n / 2; n4 = n / 4; n8 = n / 8;
14 dc = 2 * t * t; ds = sqrt(dc * (2 - dc));
15 t = 2 * dc; c = sintbl[n4] = 1; s = sintbl[0] = 0;
16 for (i = 1; i < n8; i++) {
19 sintbl[i] = s; sintbl[n4 - i] = c;
21 if (n8 != 0) sintbl[n8] = sqrt(0.5);
22 for (i = 0; i < n4; i++)
23 sintbl[n2 - i] = sintbl[i];
24 for (i = 0; i < n2 + n4; i++)
25 sintbl[i + n2] = - sintbl[i];
30 static void make_bitrev(int n, int bitrev[])
34 n2 = n / 2; i = j = 0;
39 while (k <= j) { j -= k; k /= 2; }
45 int fft(int n, Fukusosuu *x, int dum)
47 static int last_n = 0; /* {\tt n} */
48 static int *bitrev = NULL; /* */
49 static float *sintbl = NULL; /* */
50 int i, j, k, ik, h, d, k2, n4, inverse;
51 float t, s, c, dR, dI;
55 n = -n; inverse = 1; /* */
58 if (n != last_n || n == 0) {
69 if (n == 0) return 0; /* */
70 sintbl = (float *)malloc((n + n4) * sizeof(float));
71 bitrev = (int *)malloc(n * sizeof(int));
72 #else /* by T. Nishiyama */
73 sintbl = realloc(sintbl, (n + n4) * sizeof(float));
74 bitrev = realloc(bitrev, n * sizeof(int));
76 if (sintbl == NULL || bitrev == NULL) {
77 fprintf(stderr, "\n"); return 1;
79 make_sintbl(n, sintbl);
80 make_bitrev(n, bitrev);
82 for (i = 0; i < n; i++) { /* */
85 t = x[i].R; x[i].R = x[j].R; x[j].R = t;
86 t = x[i].I; x[i].I = x[j].I; x[j].I = t;
89 for (k = 1; k < n; k = k2) { /* */
91 fprintf( stderr, "%d / %d\n", k, n );
93 h = 0; k2 = k + k; d = n / k2;
94 for (j = 0; j < k; j++) {
97 fprintf( stderr, "%d / %d\r", j, k );
100 if (inverse) s = - sintbl[h];
102 for (i = j; i < n; i += k2) {
104 if( k>=4194000 ) fprintf( stderr, "in loop %d - %d < %d, k2=%d\r", j, i, n, k2 );
107 dR = s * x[ik].I + c * x[ik].R;
108 dI = c * x[ik].I - s * x[ik].R;
109 x[ik].R = x[i].R - dR; x[i].R += dR;
110 x[ik].I = x[i].I - dI; x[i].I += dI;
115 if (! inverse) /* n */
116 for (i = 0; i < n; i++) { x[i].R /= n; x[i].I /= n; }