11     data length :power of 2
13     radix       :4, 2
35             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
37             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
42             cdft(2*n, 1, a, ip, w);
45             cdft(2*n, -1, a, ip, w);
47         2*n            :data length (int)
48                         n >= 1, n = power of 2
49         a[0...2*n-1]   :input/output data (double *)
51                             a[2*j] = Re(x[j]), 
52                             a[2*j+1] = Im(x[j]), 0<=j<n
54                             a[2*k] = Re(X[k]), 
55                             a[2*k+1] = Im(X[k]), 0<=k<n
57                         length of ip >= 2+sqrt(n)
60                             2+(1<<(int)(log(n+0.5)/log(2))/2).
62         w[0...n/2-1]   :cos/sin table (double *)
66             cdft(2*n, -1, a, ip, w);
68             cdft(2*n, 1, a, ip, w);
69             for (j = 0; j <= 2 * n - 1; j++) {
78             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
79             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
81             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
82                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
83                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
93                         n >= 2, n = power of 2
97                                 a[2*k] = R[k], 0<=k<n/2
98                                 a[2*k+1] = I[k], 0<k<n/2
99                                 a[1] = R[n/2]
102                                 a[2*j] = R[j], 0<=j<n/2
103                                 a[2*j+1] = I[j], 0<j<n/2
104                                 a[1] = R[n/2]
106                         length of ip >= 2+sqrt(n/2)
109                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
111         w[0...n/2-1]   :cos/sin table (double *)
127             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
129             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
139                         n >= 2, n = power of 2
142                             a[k] = C[k], 0<=k<n
144                         length of ip >= 2+sqrt(n/2)
147                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
166             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
168             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
178                         n >= 2, n = power of 2
185                                 a[k] = S[k], 0<=k<n
188                                 a[k] = S[k], 0<k<n
191                         length of ip >= 2+sqrt(n/2)
194                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
212         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
218                         n >= 2, n = power of 2
221                             a[k] = C[k], 0<=k<=n
222         t[0...n/2]     :work area (double *)
224                         length of ip >= 2+sqrt(n/4)
227                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
248         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
254                         n >= 2, n = power of 2
257                             a[k] = S[k], 0<k<n
259         t[0...n/2-1]   :work area (double *)
261                         length of ip >= 2+sqrt(n/4)
264                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
328     if (n > (ip[0] << 2)) {  in cdft()
329         makewt(n >> 2, ip, w);  in cdft()
333             bitrv2(n, ip + 2, a);  in cdft()
336             bitrv2conj(n, ip + 2, a);  in cdft()
354     if (n > (nw << 2)) {  in rdft()
355         nw = n >> 2;  in rdft()
359     if (n > (nc << 2)) {  in rdft()
360         nc = n >> 2;  in rdft()
365             bitrv2(n, ip + 2, a);  in rdft()
379             bitrv2(n, ip + 2, a);  in rdft()
397     if (n > (nw << 2)) {  in ddct()
398         nw = n >> 2;  in ddct()
408         for (j = n - 2; j >= 2; j -= 2) {  in ddct()
416             bitrv2(n, ip + 2, a);  in ddct()
425             bitrv2(n, ip + 2, a);  in ddct()
433         for (j = 2; j < n; j += 2) {  in ddct()
451     if (n > (nw << 2)) {  in ddst()
452         nw = n >> 2;  in ddst()
462         for (j = n - 2; j >= 2; j -= 2) {  in ddst()
470             bitrv2(n, ip + 2, a);  in ddst()
479             bitrv2(n, ip + 2, a);  in ddst()
487         for (j = 2; j < n; j += 2) {  in ddst()
498     int j, k, l, m, mh, nw, nc;  in dfct()  local
520     if (n > 2) {  in dfct()
523             k = m - j;  in dfct()
526             yr = a[k] - a[n - k];  in dfct()
527             yi = a[k] + a[n - k];  in dfct()
529             a[k] = yr;  in dfct()
531             t[k] = xi + yi;  in dfct()
537             bitrv2(m, ip + 2, a);  in dfct()
545         for (j = m - 2; j >= 2; j -= 2) {  in dfct()
546             a[2 * j + 1] = a[j] + a[j + 1];  in dfct()
547             a[2 * j - 1] = a[j] - a[j + 1];  in dfct()
549         l = 2;  in dfct()
551         while (m >= 2) {  in dfct()
554                 bitrv2(m, ip + 2, t);  in dfct()
562             k = 0;  in dfct()
563             for (j = 2; j < m; j += 2) {  in dfct()
564                 k += l << 2;  in dfct()
565                 a[k - l] = t[j] - t[j + 1];  in dfct()
566                 a[k + l] = t[j] + t[j + 1];  in dfct()
571                 k = m - j;  in dfct()
572                 t[j] = t[m + k] - t[m + j];  in dfct()
573                 t[k] = t[m + k] + t[m + j];  in dfct()
579         a[n] = t[2] - t[1];  in dfct()
580         a[0] = t[2] + t[1];  in dfct()
583         a[2] = t[0];  in dfct()
591     int j, k, l, m, mh, nw, nc;  in dfst()  local
607     if (n > 2) {  in dfst()
611             k = m - j;  in dfst()
614             yr = a[k] + a[n - k];  in dfst()
615             yi = a[k] - a[n - k];  in dfst()
617             a[k] = yr;  in dfst()
619             t[k] = xi - yi;  in dfst()
626             bitrv2(m, ip + 2, a);  in dfst()
634         for (j = m - 2; j >= 2; j -= 2) {  in dfst()
635             a[2 * j + 1] = a[j] - a[j + 1];  in dfst()
636             a[2 * j - 1] = -a[j] - a[j + 1];  in dfst()
638         l = 2;  in dfst()
640         while (m >= 2) {  in dfst()
643                 bitrv2(m, ip + 2, t);  in dfst()
651             k = 0;  in dfst()
652             for (j = 2; j < m; j += 2) {  in dfst()
653                 k += l << 2;  in dfst()
654                 a[k - l] = -t[j] - t[j + 1];  in dfst()
655                 a[k + l] = t[j] - t[j + 1];  in dfst()
660                 k = m - j;  in dfst()
661                 t[j] = t[m + k] + t[m + j];  in dfst()
662                 t[k] = t[m + k] - t[m + j];  in dfst()
683     if (nw > 2) {  in makewt()
690         if (nwh > 2) {  in makewt()
691             for (j = 2; j < nwh; j += 2) {  in makewt()
699             bitrv2(nw, ip + 2, w);  in makewt()
729     int j, j1, k, k1, l, m, m2, ip[256];  in bitrv2()  local
743     m2 = 2 * m;  in bitrv2()
745         for (k = 0; k < m; k++) {  in bitrv2()
746             for (j = 0; j < k; j++) {  in bitrv2()
747                 j1 = 2 * j + ip[k];  in bitrv2()
748                 k1 = 2 * k + ip[j];  in bitrv2()
758                 k1 += 2 * m2;  in bitrv2()
778                 k1 += 2 * m2;  in bitrv2()
788             j1 = 2 * k + m2 + ip[k];  in bitrv2()
800         for (k = 1; k < m; k++) {  in bitrv2()
801             for (j = 0; j < k; j++) {  in bitrv2()
802                 j1 = 2 * j + ip[k];  in bitrv2()
803                 k1 = 2 * k + ip[j];  in bitrv2()
830     int j, j1, k, k1, l, m, m2, ip[256];  in bitrv2conj()  local
844     m2 = 2 * m;  in bitrv2conj()
846         for (k = 0; k < m; k++) {  in bitrv2conj()
847             for (j = 0; j < k; j++) {  in bitrv2conj()
848                 j1 = 2 * j + ip[k];  in bitrv2conj()
849                 k1 = 2 * k + ip[j];  in bitrv2conj()
859                 k1 += 2 * m2;  in bitrv2conj()
879                 k1 += 2 * m2;  in bitrv2conj()
889             k1 = 2 * k + ip[k];  in bitrv2conj()
907         for (k = 1; k < m; k++) {  in bitrv2conj()
908             for (j = 0; j < k; j++) {  in bitrv2conj()
909                 j1 = 2 * j + ip[k];  in bitrv2conj()
910                 k1 = 2 * k + ip[j];  in bitrv2conj()
930             k1 = 2 * k + ip[k];  in bitrv2conj()
943     l = 2;  in cftfsub()
947         while ((l << 2) < n) {  in cftfsub()
949             l <<= 2;  in cftfsub()
952     if ((l << 2) == n) {  in cftfsub()
953         for (j = 0; j < l; j += 2) {  in cftfsub()
975         for (j = 0; j < l; j += 2) {  in cftfsub()
993     l = 2;  in cftbsub()
997         while ((l << 2) < n) {  in cftbsub()
999             l <<= 2;  in cftbsub()
1002     if ((l << 2) == n) {  in cftbsub()
1003         for (j = 0; j < l; j += 2) {  in cftbsub()
1025         for (j = 0; j < l; j += 2) {  in cftbsub()
1044     x0r = a[0] + a[2];  in cft1st()
1046     x1r = a[0] - a[2];  in cft1st()
1056     a[2] = x1r - x3i;  in cft1st()
1060     wk1r = w[2];  in cft1st()
1083         k1 += 2;  in cft1st()
1084         k2 = 2 * k1;  in cft1st()
1089         wk3r = wk1r - 2 * wk2i * wk1i;  in cft1st()
1090         wk3i = 2 * wk2i * wk1r - wk1i;  in cft1st()
1091         x0r = a[j] + a[j + 2];  in cft1st()
1093         x1r = a[j] - a[j + 2];  in cft1st()
1107         a[j + 2] = wk1r * x0r - wk1i * x0i;  in cft1st()
1113         wk1r = w[k2 + 2];  in cft1st()
1115         wk3r = wk1r - 2 * wk2r * wk1i;  in cft1st()
1116         wk3i = 2 * wk2r * wk1r - wk1i;  in cft1st()
1145     int j, j1, j2, j3, k, k1, k2, m, m2;  in cftmdl()  local
1149     m = l << 2;  in cftmdl()
1150     for (j = 0; j < l; j += 2) {  in cftmdl()
1171     wk1r = w[2];  in cftmdl()
1172     for (j = m; j < l + m; j += 2) {  in cftmdl()
1198     m2 = 2 * m;  in cftmdl()
1199     for (k = m2; k < n; k += m2) {  in cftmdl()
1200         k1 += 2;  in cftmdl()
1201         k2 = 2 * k1;  in cftmdl()
1206         wk3r = wk1r - 2 * wk2i * wk1i;  in cftmdl()
1207         wk3i = 2 * wk2i * wk1r - wk1i;  in cftmdl()
1208         for (j = k; j < l + k; j += 2) {  in cftmdl()
1235         wk1r = w[k2 + 2];  in cftmdl()
1237         wk3r = wk1r - 2 * wk2r * wk1i;  in cftmdl()
1238         wk3i = 2 * wk2r * wk1r - wk1i;  in cftmdl()
1239         for (j = k + m; j < l + (k + m); j += 2) {  in cftmdl()
1272     int j, k, kk, ks, m;  in rftfsub()  local
1276     ks = 2 * nc / m;  in rftfsub()
1278     for (j = 2; j < m; j += 2) {  in rftfsub()
1279         k = n - j;  in rftfsub()
1283         xr = a[j] - a[k];  in rftfsub()
1284         xi = a[j + 1] + a[k + 1];  in rftfsub()
1289         a[k] += yr;  in rftfsub()
1290         a[k + 1] -= yi;  in rftfsub()
1297     int j, k, kk, ks, m;  in rftbsub()  local
1302     ks = 2 * nc / m;  in rftbsub()
1304     for (j = 2; j < m; j += 2) {  in rftbsub()
1305         k = n - j;  in rftbsub()
1309         xr = a[j] - a[k];  in rftbsub()
1310         xi = a[j + 1] + a[k + 1];  in rftbsub()
1315         a[k] += yr;  in rftbsub()
1316         a[k + 1] = yi - a[k + 1];  in rftbsub()
1324     int j, k, kk, ks, m;  in dctsub()  local
1331         k = n - j;  in dctsub()
1335         xr = wki * a[j] - wkr * a[k];  in dctsub()
1336         a[j] = wkr * a[j] + wki * a[k];  in dctsub()
1337         a[k] = xr;  in dctsub()
1345     int j, k, kk, ks, m;  in dstsub()  local
1352         k = n - j;  in dstsub()
1356         xr = wki * a[k] - wkr * a[j];  in dstsub()
1357         a[k] = wkr * a[k] + wki * a[j];  in dstsub()