File indexing completed on 2024-12-22 04:18:03

0001 /*
0002     Copyright(C) 1996-2001 Takuya OOURA
0003     email: ooura@mmm.t.u-tokyo.ac.jp
0004     download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
0005     You may use, copy, modify this code for any purpose and
0006     without fee.
0007 
0008 */
0009 
0010 /*
0011 Fast Fourier/Cosine/Sine Transform
0012     dimension   :one
0013     data length :power of 2
0014     decimation  :frequency
0015     radix       :split-radix
0016     data        :inplace
0017     table       :not use
0018 functions
0019     cdft: Complex Discrete Fourier Transform
0020     rdft: Real Discrete Fourier Transform
0021     ddct: Discrete Cosine Transform
0022     ddst: Discrete Sine Transform
0023     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
0024     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
0025 function prototypes
0026     void cdft(int, int, double *);
0027     void rdft(int, int, double *);
0028     void ddct(int, int, double *);
0029     void ddst(int, int, double *);
0030     void dfct(int, double *);
0031     void dfst(int, double *);
0032 macro definitions
0033     USE_CDFT_PTHREADS : default=not defined
0034         CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
0035         CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
0036     USE_CDFT_WINTHREADS : default=not defined
0037         CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
0038         CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288
0039 
0040 
0041 -------- Complex DFT (Discrete Fourier Transform) --------
0042     [definition]
0043         <case1>
0044             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
0045         <case2>
0046             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
0047         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
0048     [usage]
0049         <case1>
0050             cdft(2*n, 1, a);
0051         <case2>
0052             cdft(2*n, -1, a);
0053     [parameters]
0054         2*n            :data length (int)
0055                         n >= 1, n = power of 2
0056         a[0...2*n-1]   :input/output data (double *)
0057                         input data
0058                             a[2*j] = Re(x[j]), 
0059                             a[2*j+1] = Im(x[j]), 0<=j<n
0060                         output data
0061                             a[2*k] = Re(X[k]), 
0062                             a[2*k+1] = Im(X[k]), 0<=k<n
0063     [remark]
0064         Inverse of 
0065             cdft(2*n, -1, a);
0066         is 
0067             cdft(2*n, 1, a);
0068             for (j = 0; j <= 2 * n - 1; j++) {
0069                 a[j] *= 1.0 / n;
0070             }
0071         .
0072 
0073 
0074 -------- Real DFT / Inverse of Real DFT --------
0075     [definition]
0076         <case1> RDFT
0077             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
0078             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
0079         <case2> IRDFT (excluding scale)
0080             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
0081                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
0082                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
0083     [usage]
0084         <case1>
0085             rdft(n, 1, a);
0086         <case2>
0087             rdft(n, -1, a);
0088     [parameters]
0089         n              :data length (int)
0090                         n >= 2, n = power of 2
0091         a[0...n-1]     :input/output data (double *)
0092                         <case1>
0093                             output data
0094                                 a[2*k] = R[k], 0<=k<n/2
0095                                 a[2*k+1] = I[k], 0<k<n/2
0096                                 a[1] = R[n/2]
0097                         <case2>
0098                             input data
0099                                 a[2*j] = R[j], 0<=j<n/2
0100                                 a[2*j+1] = I[j], 0<j<n/2
0101                                 a[1] = R[n/2]
0102     [remark]
0103         Inverse of 
0104             rdft(n, 1, a);
0105         is 
0106             rdft(n, -1, a);
0107             for (j = 0; j <= n - 1; j++) {
0108                 a[j] *= 2.0 / n;
0109             }
0110         .
0111 
0112 
0113 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
0114     [definition]
0115         <case1> IDCT (excluding scale)
0116             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
0117         <case2> DCT
0118             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
0119     [usage]
0120         <case1>
0121             ddct(n, 1, a);
0122         <case2>
0123             ddct(n, -1, a);
0124     [parameters]
0125         n              :data length (int)
0126                         n >= 2, n = power of 2
0127         a[0...n-1]     :input/output data (double *)
0128                         output data
0129                             a[k] = C[k], 0<=k<n
0130     [remark]
0131         Inverse of 
0132             ddct(n, -1, a);
0133         is 
0134             a[0] *= 0.5;
0135             ddct(n, 1, a);
0136             for (j = 0; j <= n - 1; j++) {
0137                 a[j] *= 2.0 / n;
0138             }
0139         .
0140 
0141 
0142 -------- DST (Discrete Sine Transform) / Inverse of DST --------
0143     [definition]
0144         <case1> IDST (excluding scale)
0145             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
0146         <case2> DST
0147             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
0148     [usage]
0149         <case1>
0150             ddst(n, 1, a);
0151         <case2>
0152             ddst(n, -1, a);
0153     [parameters]
0154         n              :data length (int)
0155                         n >= 2, n = power of 2
0156         a[0...n-1]     :input/output data (double *)
0157                         <case1>
0158                             input data
0159                                 a[j] = A[j], 0<j<n
0160                                 a[0] = A[n]
0161                             output data
0162                                 a[k] = S[k], 0<=k<n
0163                         <case2>
0164                             output data
0165                                 a[k] = S[k], 0<k<n
0166                                 a[0] = S[n]
0167     [remark]
0168         Inverse of 
0169             ddst(n, -1, a);
0170         is 
0171             a[0] *= 0.5;
0172             ddst(n, 1, a);
0173             for (j = 0; j <= n - 1; j++) {
0174                 a[j] *= 2.0 / n;
0175             }
0176         .
0177 
0178 
0179 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
0180     [definition]
0181         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
0182     [usage]
0183         dfct(n, a);
0184     [parameters]
0185         n              :data length - 1 (int)
0186                         n >= 2, n = power of 2
0187         a[0...n]       :input/output data (double *)
0188                         output data
0189                             a[k] = C[k], 0<=k<=n
0190     [remark]
0191         Inverse of 
0192             a[0] *= 0.5;
0193             a[n] *= 0.5;
0194             dfct(n, a);
0195         is 
0196             a[0] *= 0.5;
0197             a[n] *= 0.5;
0198             dfct(n, a);
0199             for (j = 0; j <= n; j++) {
0200                 a[j] *= 2.0 / n;
0201             }
0202         .
0203 
0204 
0205 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
0206     [definition]
0207         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
0208     [usage]
0209         dfst(n, a);
0210     [parameters]
0211         n              :data length + 1 (int)
0212                         n >= 2, n = power of 2
0213         a[0...n-1]     :input/output data (double *)
0214                         output data
0215                             a[k] = S[k], 0<k<n
0216                         (a[0] is used for work area)
0217     [remark]
0218         Inverse of 
0219             dfst(n, a);
0220         is 
0221             dfst(n, a);
0222             for (j = 1; j <= n - 1; j++) {
0223                 a[j] *= 2.0 / n;
0224             }
0225         .
0226 */
0227 
0228 
0229 static void cftfsub(int n, double *a);
0230 static void cftbsub(int n, double *a);
0231 static void rftfsub(int n, double *a);
0232 static void rftbsub(int n, double *a);
0233 static void dctsub(int n, double *a);
0234 static void dctsub4(int n, double *a);
0235 static void ddct(int n, int isgn, double *a);
0236 static void bitrv1(int n, double *a);
0237 static void dstsub(int n, double *a);
0238 static void dstsub4(int n, double *a);
0239 static void bitrv2(int n, double *a);
0240 static void bitrv216(double *a);
0241 static void bitrv208(double *a);
0242 static void cftmdl1(int n, double *a);
0243 static void cftrec4(int n, double *a);
0244 static void cftleaf(int n, int isplt, double *a);
0245 static void cftfx41(int n, double *a);
0246 static void cftf161(double *a);
0247 static void cftf081(double *a);
0248 static void cftf040(double *a);
0249 static void cftx020(double *a);
0250 #ifdef USE_CDFT_THREADS
0251 static void cftrec4_th(int n, double *a);
0252 #endif /* USE_CDFT_THREADS */
0253 static void bitrv2conj(int n, double *a);
0254 static void bitrv216neg(double *a);
0255 static void bitrv208neg(double *a);
0256 static void cftb1st(int n, double *a);
0257 static void cftrec4(int n, double *a);
0258 static void cftleaf(int n, int isplt, double *a);
0259 static void cftfx41(int n, double *a);
0260 static void cftf161(double *a);
0261 static void cftf081(double *a);
0262 static void cftb040(double *a);
0263 static void cftx020(double *a);
0264 static int cfttree(int n, int j, int k, double *a);
0265 static void cftleaf(int n, int isplt, double *a);
0266 static void cftmdl1(int n, double *a);
0267 static void cftmdl2(int n, double *a);
0268 static void cftf162(double *a);
0269 static void cftf081(double *a);
0270 static void cftf082(double *a);
0271 
0272 void cdft(int n, int isgn, double *a);
0273 void rdft(int n, int isgn, double *a);
0274 void dfct(int n, double *a);
0275 void dfst(int n, double *a);
0276 
0277 void cdft(int n, int isgn, double *a)
0278 {
0279     
0280     if (isgn >= 0) {
0281         cftfsub(n, a);
0282     } else {
0283         cftbsub(n, a);
0284     }
0285 }
0286 
0287 
0288 void rdft(int n, int isgn, double *a)
0289 {
0290     double xi;
0291     
0292     if (isgn >= 0) {
0293         if (n > 4) {
0294             cftfsub(n, a);
0295             rftfsub(n, a);
0296         } else if (n == 4) {
0297             cftfsub(n, a);
0298         }
0299         xi = a[0] - a[1];
0300         a[0] += a[1];
0301         a[1] = xi;
0302     } else {
0303         a[1] = 0.5 * (a[0] - a[1]);
0304         a[0] -= a[1];
0305         if (n > 4) {
0306             rftbsub(n, a);
0307             cftbsub(n, a);
0308         } else if (n == 4) {
0309             cftbsub(n, a);
0310         }
0311     }
0312 }
0313 
0314 
0315 static void ddct(int n, int isgn, double *a)
0316 {
0317     int j;
0318     double xr;
0319     
0320     if (isgn < 0) {
0321         xr = a[n - 1];
0322         for (j = n - 2; j >= 2; j -= 2) {
0323             a[j + 1] = a[j] - a[j - 1];
0324             a[j] += a[j - 1];
0325         }
0326         a[1] = a[0] - xr;
0327         a[0] += xr;
0328         if (n > 4) {
0329             rftbsub(n, a);
0330             cftbsub(n, a);
0331         } else if (n == 4) {
0332             cftbsub(n, a);
0333         }
0334     }
0335     if (n > 4) {
0336         dctsub(n, a);
0337     } else {
0338         dctsub4(n, a);
0339     }
0340     if (isgn >= 0) {
0341         if (n > 4) {
0342             cftfsub(n, a);
0343             rftfsub(n, a);
0344         } else if (n == 4) {
0345             cftfsub(n, a);
0346         }
0347         xr = a[0] - a[1];
0348         a[0] += a[1];
0349         for (j = 2; j < n; j += 2) {
0350             a[j - 1] = a[j] - a[j + 1];
0351             a[j] += a[j + 1];
0352         }
0353         a[n - 1] = xr;
0354     }
0355 }
0356 
0357 
0358 static void ddst(int n, int isgn, double *a)
0359 {
0360     int j;
0361     double xr;
0362     
0363     if (isgn < 0) {
0364         xr = a[n - 1];
0365         for (j = n - 2; j >= 2; j -= 2) {
0366             a[j + 1] = -a[j] - a[j - 1];
0367             a[j] -= a[j - 1];
0368         }
0369         a[1] = a[0] + xr;
0370         a[0] -= xr;
0371         if (n > 4) {
0372             rftbsub(n, a);
0373             cftbsub(n, a);
0374         } else if (n == 4) {
0375             cftbsub(n, a);
0376         }
0377     }
0378     if (n > 4) {
0379         dstsub(n, a);
0380     } else {
0381         dstsub4(n, a);
0382     }
0383     if (isgn >= 0) {
0384         if (n > 4) {
0385             cftfsub(n, a);
0386             rftfsub(n, a);
0387         } else if (n == 4) {
0388             cftfsub(n, a);
0389         }
0390         xr = a[0] - a[1];
0391         a[0] += a[1];
0392         for (j = 2; j < n; j += 2) {
0393             a[j - 1] = -a[j] - a[j + 1];
0394             a[j] -= a[j + 1];
0395         }
0396         a[n - 1] = -xr;
0397     }
0398 }
0399 
0400 
0401 void dfct(int n, double *a)
0402 {
0403     int j, k, m, mh;
0404     double xr, xi, yr, yi, an;
0405     
0406     m = n >> 1;
0407     for (j = 0; j < m; j++) {
0408         k = n - j;
0409         xr = a[j] + a[k];
0410         a[j] -= a[k];
0411         a[k] = xr;
0412     }
0413     an = a[n];
0414     while (m >= 2) {
0415         ddct(m, 1, a);
0416         if (m > 2) {
0417             bitrv1(m, a);
0418         }
0419         mh = m >> 1;
0420         xi = a[m];
0421         a[m] = a[0];
0422         a[0] = an - xi;
0423         an += xi;
0424         for (j = 1; j < mh; j++) {
0425             k = m - j;
0426             xr = a[m + k];
0427             xi = a[m + j];
0428             yr = a[j];
0429             yi = a[k];
0430             a[m + j] = yr;
0431             a[m + k] = yi;
0432             a[j] = xr - xi;
0433             a[k] = xr + xi;
0434         }
0435         xr = a[mh];
0436         a[mh] = a[m + mh];
0437         a[m + mh] = xr;
0438         m = mh;
0439     }
0440     xi = a[1];
0441     a[1] = a[0];
0442     a[0] = an + xi;
0443     a[n] = an - xi;
0444     if (n > 2) {
0445         bitrv1(n, a);
0446     }
0447 }
0448 
0449 
0450 void dfst(int n, double *a)
0451 {
0452     int j, k, m, mh;
0453     double xr, xi, yr, yi;
0454     
0455     m = n >> 1;
0456     for (j = 1; j < m; j++) {
0457         k = n - j;
0458         xr = a[j] - a[k];
0459         a[j] += a[k];
0460         a[k] = xr;
0461     }
0462     a[0] = a[m];
0463     while (m >= 2) {
0464         ddst(m, 1, a);
0465         if (m > 2) {
0466             bitrv1(m, a);
0467         }
0468         mh = m >> 1;
0469         for (j = 1; j < mh; j++) {
0470             k = m - j;
0471             xr = a[m + k];
0472             xi = a[m + j];
0473             yr = a[j];
0474             yi = a[k];
0475             a[m + j] = yr;
0476             a[m + k] = yi;
0477             a[j] = xr + xi;
0478             a[k] = xr - xi;
0479         }
0480         a[m] = a[0];
0481         a[0] = a[m + mh];
0482         a[m + mh] = a[mh];
0483         m = mh;
0484     }
0485     a[1] = a[0];
0486     a[0] = 0;
0487     if (n > 2) {
0488         bitrv1(n, a);
0489     }
0490 }
0491 
0492 
0493 /* -------- child routines -------- */
0494 
0495 
0496 #include <math.h>
0497 #ifndef M_PI_2
0498 #define M_PI_2      1.570796326794896619231321691639751442098584699687
0499 #endif
0500 #ifndef WR5000  /* cos(M_PI_2*0.5000) */
0501 #define WR5000      0.707106781186547524400844362104849039284835937688
0502 #endif
0503 #ifndef WR2500  /* cos(M_PI_2*0.2500) */
0504 #define WR2500      0.923879532511286756128183189396788286822416625863
0505 #endif
0506 #ifndef WI2500  /* sin(M_PI_2*0.2500) */
0507 #define WI2500      0.382683432365089771728459984030398866761344562485
0508 #endif
0509 #ifndef WR1250  /* cos(M_PI_2*0.1250) */
0510 #define WR1250      0.980785280403230449126182236134239036973933730893
0511 #endif
0512 #ifndef WI1250  /* sin(M_PI_2*0.1250) */
0513 #define WI1250      0.195090322016128267848284868477022240927691617751
0514 #endif
0515 #ifndef WR3750  /* cos(M_PI_2*0.3750) */
0516 #define WR3750      0.831469612302545237078788377617905756738560811987
0517 #endif
0518 #ifndef WI3750  /* sin(M_PI_2*0.3750) */
0519 #define WI3750      0.555570233019602224742830813948532874374937190754
0520 #endif
0521 
0522 
0523 #ifdef USE_CDFT_PTHREADS
0524 #define USE_CDFT_THREADS
0525 #ifndef CDFT_THREADS_BEGIN_N
0526 #define CDFT_THREADS_BEGIN_N 8192
0527 #endif
0528 #ifndef CDFT_4THREADS_BEGIN_N
0529 #define CDFT_4THREADS_BEGIN_N 65536
0530 #endif
0531 #include <pthread.h>
0532 #include <stdio.h>
0533 #include <stdlib.h>
0534 #define cdft_thread_t pthread_t
0535 #define cdft_thread_create(thp,func,argp) { \
0536     if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
0537         fprintf(stderr, "cdft thread error\n"); \
0538         exit(1); \
0539     } \
0540 }
0541 #define cdft_thread_wait(th) { \
0542     if (pthread_join(th, NULL) != 0) { \
0543         fprintf(stderr, "cdft thread error\n"); \
0544         exit(1); \
0545     } \
0546 }
0547 #endif /* USE_CDFT_PTHREADS */
0548 
0549 
0550 #ifdef USE_CDFT_WINTHREADS
0551 #define USE_CDFT_THREADS
0552 #ifndef CDFT_THREADS_BEGIN_N
0553 #define CDFT_THREADS_BEGIN_N 32768
0554 #endif
0555 #ifndef CDFT_4THREADS_BEGIN_N
0556 #define CDFT_4THREADS_BEGIN_N 524288
0557 #endif
0558 #include <windows.h>
0559 #include <stdio.h>
0560 #include <stdlib.h>
0561 #define cdft_thread_t HANDLE
0562 #define cdft_thread_create(thp,func,argp) { \
0563     DWORD thid; \
0564     *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
0565     if (*(thp) == 0) { \
0566         fprintf(stderr, "cdft thread error\n"); \
0567         exit(1); \
0568     } \
0569 }
0570 #define cdft_thread_wait(th) { \
0571     WaitForSingleObject(th, INFINITE); \
0572     CloseHandle(th); \
0573 }
0574 #endif /* USE_CDFT_WINTHREADS */
0575 
0576 
0577 #ifndef CDFT_LOOP_DIV  /* control of the CDFT's speed & tolerance */
0578 #define CDFT_LOOP_DIV 32
0579 #endif
0580 
0581 #ifndef RDFT_LOOP_DIV  /* control of the RDFT's speed & tolerance */
0582 #define RDFT_LOOP_DIV 64
0583 #endif
0584 
0585 #ifndef DCST_LOOP_DIV  /* control of the DCT,DST's speed & tolerance */
0586 #define DCST_LOOP_DIV 64
0587 #endif
0588 
0589 
0590 static void cftfsub(int n, double *a)
0591 {
0592     
0593     if (n > 8) {
0594         if (n > 32) {
0595             cftmdl1(n, a);
0596 #ifdef USE_CDFT_THREADS
0597             if (n > CDFT_THREADS_BEGIN_N) {
0598                 cftrec4_th(n, a);
0599             } else 
0600 #endif /* USE_CDFT_THREADS */
0601             if (n > 512) {
0602                 cftrec4(n, a);
0603             } else if (n > 128) {
0604                 cftleaf(n, 1, a);
0605             } else {
0606                 cftfx41(n, a);
0607             }
0608             bitrv2(n, a);
0609         } else if (n == 32) {
0610             cftf161(a);
0611             bitrv216(a);
0612         } else {
0613             cftf081(a);
0614             bitrv208(a);
0615         }
0616     } else if (n == 8) {
0617         cftf040(a);
0618     } else if (n == 4) {
0619         cftx020(a);
0620     }
0621 }
0622 
0623 
0624 static void cftbsub(int n, double *a)
0625 {
0626     if (n > 8) {
0627         if (n > 32) {
0628             cftb1st(n, a);
0629 #ifdef USE_CDFT_THREADS
0630             if (n > CDFT_THREADS_BEGIN_N) {
0631                 cftrec4_th(n, a);
0632             } else 
0633 #endif /* USE_CDFT_THREADS */
0634             if (n > 512) {
0635                 cftrec4(n, a);
0636             } else if (n > 128) {
0637                 cftleaf(n, 1, a);
0638             } else {
0639                 cftfx41(n, a);
0640             }
0641             bitrv2conj(n, a);
0642         } else if (n == 32) {
0643             cftf161(a);
0644             bitrv216neg(a);
0645         } else {
0646             cftf081(a);
0647             bitrv208neg(a);
0648         }
0649     } else if (n == 8) {
0650         cftb040(a);
0651     } else if (n == 4) {
0652         cftx020(a);
0653     }
0654 }
0655 
0656 
0657 void bitrv2(int n, double *a)
0658 {
0659     int _j0, k0, _j1, k1, l, m, i, j, k, nh;
0660     double xr, xi, yr, yi;
0661     
0662     m = 4;
0663     for (l = n >> 2; l > 8; l >>= 2) {
0664         m <<= 1;
0665     }
0666     nh = n >> 1;
0667     if (l == 8) {
0668         _j0 = 0;
0669         for (k0 = 0; k0 < m; k0 += 4) {
0670             k = k0;
0671             for (j = _j0; j < _j0 + k0; j += 4) {
0672                 xr = a[j];
0673                 xi = a[j + 1];
0674                 yr = a[k];
0675                 yi = a[k + 1];
0676                 a[j] = yr;
0677                 a[j + 1] = yi;
0678                 a[k] = xr;
0679                 a[k + 1] = xi;
0680                 _j1 = j + m;
0681                 k1 = k + 2 * m;
0682                 xr = a[_j1];
0683                 xi = a[_j1 + 1];
0684                 yr = a[k1];
0685                 yi = a[k1 + 1];
0686                 a[_j1] = yr;
0687                 a[_j1 + 1] = yi;
0688                 a[k1] = xr;
0689                 a[k1 + 1] = xi;
0690                 _j1 += m;
0691                 k1 -= m;
0692                 xr = a[_j1];
0693                 xi = a[_j1 + 1];
0694                 yr = a[k1];
0695                 yi = a[k1 + 1];
0696                 a[_j1] = yr;
0697                 a[_j1 + 1] = yi;
0698                 a[k1] = xr;
0699                 a[k1 + 1] = xi;
0700                 _j1 += m;
0701                 k1 += 2 * m;
0702                 xr = a[_j1];
0703                 xi = a[_j1 + 1];
0704                 yr = a[k1];
0705                 yi = a[k1 + 1];
0706                 a[_j1] = yr;
0707                 a[_j1 + 1] = yi;
0708                 a[k1] = xr;
0709                 a[k1 + 1] = xi;
0710                 _j1 += nh;
0711                 k1 += 2;
0712                 xr = a[_j1];
0713                 xi = a[_j1 + 1];
0714                 yr = a[k1];
0715                 yi = a[k1 + 1];
0716                 a[_j1] = yr;
0717                 a[_j1 + 1] = yi;
0718                 a[k1] = xr;
0719                 a[k1 + 1] = xi;
0720                 _j1 -= m;
0721                 k1 -= 2 * m;
0722                 xr = a[_j1];
0723                 xi = a[_j1 + 1];
0724                 yr = a[k1];
0725                 yi = a[k1 + 1];
0726                 a[_j1] = yr;
0727                 a[_j1 + 1] = yi;
0728                 a[k1] = xr;
0729                 a[k1 + 1] = xi;
0730                 _j1 -= m;
0731                 k1 += m;
0732                 xr = a[_j1];
0733                 xi = a[_j1 + 1];
0734                 yr = a[k1];
0735                 yi = a[k1 + 1];
0736                 a[_j1] = yr;
0737                 a[_j1 + 1] = yi;
0738                 a[k1] = xr;
0739                 a[k1 + 1] = xi;
0740                 _j1 -= m;
0741                 k1 -= 2 * m;
0742                 xr = a[_j1];
0743                 xi = a[_j1 + 1];
0744                 yr = a[k1];
0745                 yi = a[k1 + 1];
0746                 a[_j1] = yr;
0747                 a[_j1 + 1] = yi;
0748                 a[k1] = xr;
0749                 a[k1 + 1] = xi;
0750                 _j1 += 2;
0751                 k1 += nh;
0752                 xr = a[_j1];
0753                 xi = a[_j1 + 1];
0754                 yr = a[k1];
0755                 yi = a[k1 + 1];
0756                 a[_j1] = yr;
0757                 a[_j1 + 1] = yi;
0758                 a[k1] = xr;
0759                 a[k1 + 1] = xi;
0760                 _j1 += m;
0761                 k1 += 2 * m;
0762                 xr = a[_j1];
0763                 xi = a[_j1 + 1];
0764                 yr = a[k1];
0765                 yi = a[k1 + 1];
0766                 a[_j1] = yr;
0767                 a[_j1 + 1] = yi;
0768                 a[k1] = xr;
0769                 a[k1 + 1] = xi;
0770                 _j1 += m;
0771                 k1 -= m;
0772                 xr = a[_j1];
0773                 xi = a[_j1 + 1];
0774                 yr = a[k1];
0775                 yi = a[k1 + 1];
0776                 a[_j1] = yr;
0777                 a[_j1 + 1] = yi;
0778                 a[k1] = xr;
0779                 a[k1 + 1] = xi;
0780                 _j1 += m;
0781                 k1 += 2 * m;
0782                 xr = a[_j1];
0783                 xi = a[_j1 + 1];
0784                 yr = a[k1];
0785                 yi = a[k1 + 1];
0786                 a[_j1] = yr;
0787                 a[_j1 + 1] = yi;
0788                 a[k1] = xr;
0789                 a[k1 + 1] = xi;
0790                 _j1 -= nh;
0791                 k1 -= 2;
0792                 xr = a[_j1];
0793                 xi = a[_j1 + 1];
0794                 yr = a[k1];
0795                 yi = a[k1 + 1];
0796                 a[_j1] = yr;
0797                 a[_j1 + 1] = yi;
0798                 a[k1] = xr;
0799                 a[k1 + 1] = xi;
0800                 _j1 -= m;
0801                 k1 -= 2 * m;
0802                 xr = a[_j1];
0803                 xi = a[_j1 + 1];
0804                 yr = a[k1];
0805                 yi = a[k1 + 1];
0806                 a[_j1] = yr;
0807                 a[_j1 + 1] = yi;
0808                 a[k1] = xr;
0809                 a[k1 + 1] = xi;
0810                 _j1 -= m;
0811                 k1 += m;
0812                 xr = a[_j1];
0813                 xi = a[_j1 + 1];
0814                 yr = a[k1];
0815                 yi = a[k1 + 1];
0816                 a[_j1] = yr;
0817                 a[_j1 + 1] = yi;
0818                 a[k1] = xr;
0819                 a[k1 + 1] = xi;
0820                 _j1 -= m;
0821                 k1 -= 2 * m;
0822                 xr = a[_j1];
0823                 xi = a[_j1 + 1];
0824                 yr = a[k1];
0825                 yi = a[k1 + 1];
0826                 a[_j1] = yr;
0827                 a[_j1 + 1] = yi;
0828                 a[k1] = xr;
0829                 a[k1 + 1] = xi;
0830                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
0831             }
0832             k1 = _j0 + k0;
0833             _j1 = k1 + 2;
0834             k1 += nh;
0835             xr = a[_j1];
0836             xi = a[_j1 + 1];
0837             yr = a[k1];
0838             yi = a[k1 + 1];
0839             a[_j1] = yr;
0840             a[_j1 + 1] = yi;
0841             a[k1] = xr;
0842             a[k1 + 1] = xi;
0843             _j1 += m;
0844             k1 += 2 * m;
0845             xr = a[_j1];
0846             xi = a[_j1 + 1];
0847             yr = a[k1];
0848             yi = a[k1 + 1];
0849             a[_j1] = yr;
0850             a[_j1 + 1] = yi;
0851             a[k1] = xr;
0852             a[k1 + 1] = xi;
0853             _j1 += m;
0854             k1 -= m;
0855             xr = a[_j1];
0856             xi = a[_j1 + 1];
0857             yr = a[k1];
0858             yi = a[k1 + 1];
0859             a[_j1] = yr;
0860             a[_j1 + 1] = yi;
0861             a[k1] = xr;
0862             a[k1 + 1] = xi;
0863             _j1 -= 2;
0864             k1 -= nh;
0865             xr = a[_j1];
0866             xi = a[_j1 + 1];
0867             yr = a[k1];
0868             yi = a[k1 + 1];
0869             a[_j1] = yr;
0870             a[_j1 + 1] = yi;
0871             a[k1] = xr;
0872             a[k1 + 1] = xi;
0873             _j1 += nh + 2;
0874             k1 += nh + 2;
0875             xr = a[_j1];
0876             xi = a[_j1 + 1];
0877             yr = a[k1];
0878             yi = a[k1 + 1];
0879             a[_j1] = yr;
0880             a[_j1 + 1] = yi;
0881             a[k1] = xr;
0882             a[k1 + 1] = xi;
0883             _j1 -= nh - m;
0884             k1 += 2 * m - 2;
0885             xr = a[_j1];
0886             xi = a[_j1 + 1];
0887             yr = a[k1];
0888             yi = a[k1 + 1];
0889             a[_j1] = yr;
0890             a[_j1 + 1] = yi;
0891             a[k1] = xr;
0892             a[k1 + 1] = xi;
0893             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
0894         }
0895     } else {
0896         _j0 = 0;
0897         for (k0 = 0; k0 < m; k0 += 4) {
0898             k = k0;
0899             for (j = _j0; j < _j0 + k0; j += 4) {
0900                 xr = a[j];
0901                 xi = a[j + 1];
0902                 yr = a[k];
0903                 yi = a[k + 1];
0904                 a[j] = yr;
0905                 a[j + 1] = yi;
0906                 a[k] = xr;
0907                 a[k + 1] = xi;
0908                 _j1 = j + m;
0909                 k1 = k + m;
0910                 xr = a[_j1];
0911                 xi = a[_j1 + 1];
0912                 yr = a[k1];
0913                 yi = a[k1 + 1];
0914                 a[_j1] = yr;
0915                 a[_j1 + 1] = yi;
0916                 a[k1] = xr;
0917                 a[k1 + 1] = xi;
0918                 _j1 += nh;
0919                 k1 += 2;
0920                 xr = a[_j1];
0921                 xi = a[_j1 + 1];
0922                 yr = a[k1];
0923                 yi = a[k1 + 1];
0924                 a[_j1] = yr;
0925                 a[_j1 + 1] = yi;
0926                 a[k1] = xr;
0927                 a[k1 + 1] = xi;
0928                 _j1 -= m;
0929                 k1 -= m;
0930                 xr = a[_j1];
0931                 xi = a[_j1 + 1];
0932                 yr = a[k1];
0933                 yi = a[k1 + 1];
0934                 a[_j1] = yr;
0935                 a[_j1 + 1] = yi;
0936                 a[k1] = xr;
0937                 a[k1 + 1] = xi;
0938                 _j1 += 2;
0939                 k1 += nh;
0940                 xr = a[_j1];
0941                 xi = a[_j1 + 1];
0942                 yr = a[k1];
0943                 yi = a[k1 + 1];
0944                 a[_j1] = yr;
0945                 a[_j1 + 1] = yi;
0946                 a[k1] = xr;
0947                 a[k1 + 1] = xi;
0948                 _j1 += m;
0949                 k1 += m;
0950                 xr = a[_j1];
0951                 xi = a[_j1 + 1];
0952                 yr = a[k1];
0953                 yi = a[k1 + 1];
0954                 a[_j1] = yr;
0955                 a[_j1 + 1] = yi;
0956                 a[k1] = xr;
0957                 a[k1 + 1] = xi;
0958                 _j1 -= nh;
0959                 k1 -= 2;
0960                 xr = a[_j1];
0961                 xi = a[_j1 + 1];
0962                 yr = a[k1];
0963                 yi = a[k1 + 1];
0964                 a[_j1] = yr;
0965                 a[_j1 + 1] = yi;
0966                 a[k1] = xr;
0967                 a[k1 + 1] = xi;
0968                 _j1 -= m;
0969                 k1 -= m;
0970                 xr = a[_j1];
0971                 xi = a[_j1 + 1];
0972                 yr = a[k1];
0973                 yi = a[k1 + 1];
0974                 a[_j1] = yr;
0975                 a[_j1 + 1] = yi;
0976                 a[k1] = xr;
0977                 a[k1 + 1] = xi;
0978                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
0979             }
0980             k1 = _j0 + k0;
0981             _j1 = k1 + 2;
0982             k1 += nh;
0983             xr = a[_j1];
0984             xi = a[_j1 + 1];
0985             yr = a[k1];
0986             yi = a[k1 + 1];
0987             a[_j1] = yr;
0988             a[_j1 + 1] = yi;
0989             a[k1] = xr;
0990             a[k1 + 1] = xi;
0991             _j1 += m;
0992             k1 += m;
0993             xr = a[_j1];
0994             xi = a[_j1 + 1];
0995             yr = a[k1];
0996             yi = a[k1 + 1];
0997             a[_j1] = yr;
0998             a[_j1 + 1] = yi;
0999             a[k1] = xr;
1000             a[k1 + 1] = xi;
1001             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
1002         }
1003     }
1004 }
1005 
1006 
1007 void bitrv2conj(int n, double *a)
1008 {
1009     int _j0, k0, _j1, k1, l, m, i, j, k, nh;
1010     double xr, xi, yr, yi;
1011     
1012     m = 4;
1013     for (l = n >> 2; l > 8; l >>= 2) {
1014         m <<= 1;
1015     }
1016     nh = n >> 1;
1017     if (l == 8) {
1018         _j0 = 0;
1019         for (k0 = 0; k0 < m; k0 += 4) {
1020             k = k0;
1021             for (j = _j0; j < _j0 + k0; j += 4) {
1022                 xr = a[j];
1023                 xi = -a[j + 1];
1024                 yr = a[k];
1025                 yi = -a[k + 1];
1026                 a[j] = yr;
1027                 a[j + 1] = yi;
1028                 a[k] = xr;
1029                 a[k + 1] = xi;
1030                 _j1 = j + m;
1031                 k1 = k + 2 * m;
1032                 xr = a[_j1];
1033                 xi = -a[_j1 + 1];
1034                 yr = a[k1];
1035                 yi = -a[k1 + 1];
1036                 a[_j1] = yr;
1037                 a[_j1 + 1] = yi;
1038                 a[k1] = xr;
1039                 a[k1 + 1] = xi;
1040                 _j1 += m;
1041                 k1 -= m;
1042                 xr = a[_j1];
1043                 xi = -a[_j1 + 1];
1044                 yr = a[k1];
1045                 yi = -a[k1 + 1];
1046                 a[_j1] = yr;
1047                 a[_j1 + 1] = yi;
1048                 a[k1] = xr;
1049                 a[k1 + 1] = xi;
1050                 _j1 += m;
1051                 k1 += 2 * m;
1052                 xr = a[_j1];
1053                 xi = -a[_j1 + 1];
1054                 yr = a[k1];
1055                 yi = -a[k1 + 1];
1056                 a[_j1] = yr;
1057                 a[_j1 + 1] = yi;
1058                 a[k1] = xr;
1059                 a[k1 + 1] = xi;
1060                 _j1 += nh;
1061                 k1 += 2;
1062                 xr = a[_j1];
1063                 xi = -a[_j1 + 1];
1064                 yr = a[k1];
1065                 yi = -a[k1 + 1];
1066                 a[_j1] = yr;
1067                 a[_j1 + 1] = yi;
1068                 a[k1] = xr;
1069                 a[k1 + 1] = xi;
1070                 _j1 -= m;
1071                 k1 -= 2 * m;
1072                 xr = a[_j1];
1073                 xi = -a[_j1 + 1];
1074                 yr = a[k1];
1075                 yi = -a[k1 + 1];
1076                 a[_j1] = yr;
1077                 a[_j1 + 1] = yi;
1078                 a[k1] = xr;
1079                 a[k1 + 1] = xi;
1080                 _j1 -= m;
1081                 k1 += m;
1082                 xr = a[_j1];
1083                 xi = -a[_j1 + 1];
1084                 yr = a[k1];
1085                 yi = -a[k1 + 1];
1086                 a[_j1] = yr;
1087                 a[_j1 + 1] = yi;
1088                 a[k1] = xr;
1089                 a[k1 + 1] = xi;
1090                 _j1 -= m;
1091                 k1 -= 2 * m;
1092                 xr = a[_j1];
1093                 xi = -a[_j1 + 1];
1094                 yr = a[k1];
1095                 yi = -a[k1 + 1];
1096                 a[_j1] = yr;
1097                 a[_j1 + 1] = yi;
1098                 a[k1] = xr;
1099                 a[k1 + 1] = xi;
1100                 _j1 += 2;
1101                 k1 += nh;
1102                 xr = a[_j1];
1103                 xi = -a[_j1 + 1];
1104                 yr = a[k1];
1105                 yi = -a[k1 + 1];
1106                 a[_j1] = yr;
1107                 a[_j1 + 1] = yi;
1108                 a[k1] = xr;
1109                 a[k1 + 1] = xi;
1110                 _j1 += m;
1111                 k1 += 2 * m;
1112                 xr = a[_j1];
1113                 xi = -a[_j1 + 1];
1114                 yr = a[k1];
1115                 yi = -a[k1 + 1];
1116                 a[_j1] = yr;
1117                 a[_j1 + 1] = yi;
1118                 a[k1] = xr;
1119                 a[k1 + 1] = xi;
1120                 _j1 += m;
1121                 k1 -= m;
1122                 xr = a[_j1];
1123                 xi = -a[_j1 + 1];
1124                 yr = a[k1];
1125                 yi = -a[k1 + 1];
1126                 a[_j1] = yr;
1127                 a[_j1 + 1] = yi;
1128                 a[k1] = xr;
1129                 a[k1 + 1] = xi;
1130                 _j1 += m;
1131                 k1 += 2 * m;
1132                 xr = a[_j1];
1133                 xi = -a[_j1 + 1];
1134                 yr = a[k1];
1135                 yi = -a[k1 + 1];
1136                 a[_j1] = yr;
1137                 a[_j1 + 1] = yi;
1138                 a[k1] = xr;
1139                 a[k1 + 1] = xi;
1140                 _j1 -= nh;
1141                 k1 -= 2;
1142                 xr = a[_j1];
1143                 xi = -a[_j1 + 1];
1144                 yr = a[k1];
1145                 yi = -a[k1 + 1];
1146                 a[_j1] = yr;
1147                 a[_j1 + 1] = yi;
1148                 a[k1] = xr;
1149                 a[k1 + 1] = xi;
1150                 _j1 -= m;
1151                 k1 -= 2 * m;
1152                 xr = a[_j1];
1153                 xi = -a[_j1 + 1];
1154                 yr = a[k1];
1155                 yi = -a[k1 + 1];
1156                 a[_j1] = yr;
1157                 a[_j1 + 1] = yi;
1158                 a[k1] = xr;
1159                 a[k1 + 1] = xi;
1160                 _j1 -= m;
1161                 k1 += m;
1162                 xr = a[_j1];
1163                 xi = -a[_j1 + 1];
1164                 yr = a[k1];
1165                 yi = -a[k1 + 1];
1166                 a[_j1] = yr;
1167                 a[_j1 + 1] = yi;
1168                 a[k1] = xr;
1169                 a[k1 + 1] = xi;
1170                 _j1 -= m;
1171                 k1 -= 2 * m;
1172                 xr = a[_j1];
1173                 xi = -a[_j1 + 1];
1174                 yr = a[k1];
1175                 yi = -a[k1 + 1];
1176                 a[_j1] = yr;
1177                 a[_j1 + 1] = yi;
1178                 a[k1] = xr;
1179                 a[k1 + 1] = xi;
1180                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
1181             }
1182             k1 = _j0 + k0;
1183             _j1 = k1 + 2;
1184             k1 += nh;
1185             a[_j1 - 1] = -a[_j1 - 1];
1186             xr = a[_j1];
1187             xi = -a[_j1 + 1];
1188             yr = a[k1];
1189             yi = -a[k1 + 1];
1190             a[_j1] = yr;
1191             a[_j1 + 1] = yi;
1192             a[k1] = xr;
1193             a[k1 + 1] = xi;
1194             a[k1 + 3] = -a[k1 + 3];
1195             _j1 += m;
1196             k1 += 2 * m;
1197             xr = a[_j1];
1198             xi = -a[_j1 + 1];
1199             yr = a[k1];
1200             yi = -a[k1 + 1];
1201             a[_j1] = yr;
1202             a[_j1 + 1] = yi;
1203             a[k1] = xr;
1204             a[k1 + 1] = xi;
1205             _j1 += m;
1206             k1 -= m;
1207             xr = a[_j1];
1208             xi = -a[_j1 + 1];
1209             yr = a[k1];
1210             yi = -a[k1 + 1];
1211             a[_j1] = yr;
1212             a[_j1 + 1] = yi;
1213             a[k1] = xr;
1214             a[k1 + 1] = xi;
1215             _j1 -= 2;
1216             k1 -= nh;
1217             xr = a[_j1];
1218             xi = -a[_j1 + 1];
1219             yr = a[k1];
1220             yi = -a[k1 + 1];
1221             a[_j1] = yr;
1222             a[_j1 + 1] = yi;
1223             a[k1] = xr;
1224             a[k1 + 1] = xi;
1225             _j1 += nh + 2;
1226             k1 += nh + 2;
1227             xr = a[_j1];
1228             xi = -a[_j1 + 1];
1229             yr = a[k1];
1230             yi = -a[k1 + 1];
1231             a[_j1] = yr;
1232             a[_j1 + 1] = yi;
1233             a[k1] = xr;
1234             a[k1 + 1] = xi;
1235             _j1 -= nh - m;
1236             k1 += 2 * m - 2;
1237             a[_j1 - 1] = -a[_j1 - 1];
1238             xr = a[_j1];
1239             xi = -a[_j1 + 1];
1240             yr = a[k1];
1241             yi = -a[k1 + 1];
1242             a[_j1] = yr;
1243             a[_j1 + 1] = yi;
1244             a[k1] = xr;
1245             a[k1 + 1] = xi;
1246             a[k1 + 3] = -a[k1 + 3];
1247             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
1248         }
1249     } else {
1250         _j0 = 0;
1251         for (k0 = 0; k0 < m; k0 += 4) {
1252             k = k0;
1253             for (j = _j0; j < _j0 + k0; j += 4) {
1254                 xr = a[j];
1255                 xi = -a[j + 1];
1256                 yr = a[k];
1257                 yi = -a[k + 1];
1258                 a[j] = yr;
1259                 a[j + 1] = yi;
1260                 a[k] = xr;
1261                 a[k + 1] = xi;
1262                 _j1 = j + m;
1263                 k1 = k + m;
1264                 xr = a[_j1];
1265                 xi = -a[_j1 + 1];
1266                 yr = a[k1];
1267                 yi = -a[k1 + 1];
1268                 a[_j1] = yr;
1269                 a[_j1 + 1] = yi;
1270                 a[k1] = xr;
1271                 a[k1 + 1] = xi;
1272                 _j1 += nh;
1273                 k1 += 2;
1274                 xr = a[_j1];
1275                 xi = -a[_j1 + 1];
1276                 yr = a[k1];
1277                 yi = -a[k1 + 1];
1278                 a[_j1] = yr;
1279                 a[_j1 + 1] = yi;
1280                 a[k1] = xr;
1281                 a[k1 + 1] = xi;
1282                 _j1 -= m;
1283                 k1 -= m;
1284                 xr = a[_j1];
1285                 xi = -a[_j1 + 1];
1286                 yr = a[k1];
1287                 yi = -a[k1 + 1];
1288                 a[_j1] = yr;
1289                 a[_j1 + 1] = yi;
1290                 a[k1] = xr;
1291                 a[k1 + 1] = xi;
1292                 _j1 += 2;
1293                 k1 += nh;
1294                 xr = a[_j1];
1295                 xi = -a[_j1 + 1];
1296                 yr = a[k1];
1297                 yi = -a[k1 + 1];
1298                 a[_j1] = yr;
1299                 a[_j1 + 1] = yi;
1300                 a[k1] = xr;
1301                 a[k1 + 1] = xi;
1302                 _j1 += m;
1303                 k1 += m;
1304                 xr = a[_j1];
1305                 xi = -a[_j1 + 1];
1306                 yr = a[k1];
1307                 yi = -a[k1 + 1];
1308                 a[_j1] = yr;
1309                 a[_j1 + 1] = yi;
1310                 a[k1] = xr;
1311                 a[k1 + 1] = xi;
1312                 _j1 -= nh;
1313                 k1 -= 2;
1314                 xr = a[_j1];
1315                 xi = -a[_j1 + 1];
1316                 yr = a[k1];
1317                 yi = -a[k1 + 1];
1318                 a[_j1] = yr;
1319                 a[_j1 + 1] = yi;
1320                 a[k1] = xr;
1321                 a[k1 + 1] = xi;
1322                 _j1 -= m;
1323                 k1 -= m;
1324                 xr = a[_j1];
1325                 xi = -a[_j1 + 1];
1326                 yr = a[k1];
1327                 yi = -a[k1 + 1];
1328                 a[_j1] = yr;
1329                 a[_j1 + 1] = yi;
1330                 a[k1] = xr;
1331                 a[k1 + 1] = xi;
1332                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
1333             }
1334             k1 = _j0 + k0;
1335             _j1 = k1 + 2;
1336             k1 += nh;
1337             a[_j1 - 1] = -a[_j1 - 1];
1338             xr = a[_j1];
1339             xi = -a[_j1 + 1];
1340             yr = a[k1];
1341             yi = -a[k1 + 1];
1342             a[_j1] = yr;
1343             a[_j1 + 1] = yi;
1344             a[k1] = xr;
1345             a[k1 + 1] = xi;
1346             a[k1 + 3] = -a[k1 + 3];
1347             _j1 += m;
1348             k1 += m;
1349             a[_j1 - 1] = -a[_j1 - 1];
1350             xr = a[_j1];
1351             xi = -a[_j1 + 1];
1352             yr = a[k1];
1353             yi = -a[k1 + 1];
1354             a[_j1] = yr;
1355             a[_j1 + 1] = yi;
1356             a[k1] = xr;
1357             a[k1 + 1] = xi;
1358             a[k1 + 3] = -a[k1 + 3];
1359             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
1360         }
1361     }
1362 }
1363 
1364 
1365 void bitrv216(double *a)
1366 {
1367     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
1368         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
1369         x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1370     
1371     x1r = a[2];
1372     x1i = a[3];
1373     x2r = a[4];
1374     x2i = a[5];
1375     x3r = a[6];
1376     x3i = a[7];
1377     x4r = a[8];
1378     x4i = a[9];
1379     x5r = a[10];
1380     x5i = a[11];
1381     x7r = a[14];
1382     x7i = a[15];
1383     x8r = a[16];
1384     x8i = a[17];
1385     x10r = a[20];
1386     x10i = a[21];
1387     x11r = a[22];
1388     x11i = a[23];
1389     x12r = a[24];
1390     x12i = a[25];
1391     x13r = a[26];
1392     x13i = a[27];
1393     x14r = a[28];
1394     x14i = a[29];
1395     a[2] = x8r;
1396     a[3] = x8i;
1397     a[4] = x4r;
1398     a[5] = x4i;
1399     a[6] = x12r;
1400     a[7] = x12i;
1401     a[8] = x2r;
1402     a[9] = x2i;
1403     a[10] = x10r;
1404     a[11] = x10i;
1405     a[14] = x14r;
1406     a[15] = x14i;
1407     a[16] = x1r;
1408     a[17] = x1i;
1409     a[20] = x5r;
1410     a[21] = x5i;
1411     a[22] = x13r;
1412     a[23] = x13i;
1413     a[24] = x3r;
1414     a[25] = x3i;
1415     a[26] = x11r;
1416     a[27] = x11i;
1417     a[28] = x7r;
1418     a[29] = x7i;
1419 }
1420 
1421 
1422 void bitrv216neg(double *a)
1423 {
1424     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
1425         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
1426         x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, 
1427         x13r, x13i, x14r, x14i, x15r, x15i;
1428     
1429     x1r = a[2];
1430     x1i = a[3];
1431     x2r = a[4];
1432     x2i = a[5];
1433     x3r = a[6];
1434     x3i = a[7];
1435     x4r = a[8];
1436     x4i = a[9];
1437     x5r = a[10];
1438     x5i = a[11];
1439     x6r = a[12];
1440     x6i = a[13];
1441     x7r = a[14];
1442     x7i = a[15];
1443     x8r = a[16];
1444     x8i = a[17];
1445     x9r = a[18];
1446     x9i = a[19];
1447     x10r = a[20];
1448     x10i = a[21];
1449     x11r = a[22];
1450     x11i = a[23];
1451     x12r = a[24];
1452     x12i = a[25];
1453     x13r = a[26];
1454     x13i = a[27];
1455     x14r = a[28];
1456     x14i = a[29];
1457     x15r = a[30];
1458     x15i = a[31];
1459     a[2] = x15r;
1460     a[3] = x15i;
1461     a[4] = x7r;
1462     a[5] = x7i;
1463     a[6] = x11r;
1464     a[7] = x11i;
1465     a[8] = x3r;
1466     a[9] = x3i;
1467     a[10] = x13r;
1468     a[11] = x13i;
1469     a[12] = x5r;
1470     a[13] = x5i;
1471     a[14] = x9r;
1472     a[15] = x9i;
1473     a[16] = x1r;
1474     a[17] = x1i;
1475     a[18] = x14r;
1476     a[19] = x14i;
1477     a[20] = x6r;
1478     a[21] = x6i;
1479     a[22] = x10r;
1480     a[23] = x10i;
1481     a[24] = x2r;
1482     a[25] = x2i;
1483     a[26] = x12r;
1484     a[27] = x12i;
1485     a[28] = x4r;
1486     a[29] = x4i;
1487     a[30] = x8r;
1488     a[31] = x8i;
1489 }
1490 
1491 
1492 void bitrv208(double *a)
1493 {
1494     double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1495     
1496     x1r = a[2];
1497     x1i = a[3];
1498     x3r = a[6];
1499     x3i = a[7];
1500     x4r = a[8];
1501     x4i = a[9];
1502     x6r = a[12];
1503     x6i = a[13];
1504     a[2] = x4r;
1505     a[3] = x4i;
1506     a[6] = x6r;
1507     a[7] = x6i;
1508     a[8] = x1r;
1509     a[9] = x1i;
1510     a[12] = x3r;
1511     a[13] = x3i;
1512 }
1513 
1514 
1515 void bitrv208neg(double *a)
1516 {
1517     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
1518         x5r, x5i, x6r, x6i, x7r, x7i;
1519     
1520     x1r = a[2];
1521     x1i = a[3];
1522     x2r = a[4];
1523     x2i = a[5];
1524     x3r = a[6];
1525     x3i = a[7];
1526     x4r = a[8];
1527     x4i = a[9];
1528     x5r = a[10];
1529     x5i = a[11];
1530     x6r = a[12];
1531     x6i = a[13];
1532     x7r = a[14];
1533     x7i = a[15];
1534     a[2] = x7r;
1535     a[3] = x7i;
1536     a[4] = x3r;
1537     a[5] = x3i;
1538     a[6] = x5r;
1539     a[7] = x5i;
1540     a[8] = x1r;
1541     a[9] = x1i;
1542     a[10] = x6r;
1543     a[11] = x6i;
1544     a[12] = x2r;
1545     a[13] = x2i;
1546     a[14] = x4r;
1547     a[15] = x4i;
1548 }
1549 
1550 
1551 void bitrv1(int n, double *a)
1552 {
1553     int _j0, k0, _j1, k1, l, m, i, j, k, nh;
1554     double x;
1555     
1556     nh = n >> 1;
1557     x = a[1];
1558     a[1] = a[nh];
1559     a[nh] = x;
1560     m = 2;
1561     for (l = n >> 2; l > 2; l >>= 2) {
1562         m <<= 1;
1563     }
1564     if (l == 2) {
1565         _j1 = m + 1;
1566         k1 = m + nh;
1567         x = a[_j1];
1568         a[_j1] = a[k1];
1569         a[k1] = x;
1570         _j0 = 0;
1571         for (k0 = 2; k0 < m; k0 += 2) {
1572             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
1573             k = k0;
1574             for (j = _j0; j < _j0 + k0; j += 2) {
1575                 x = a[j];
1576                 a[j] = a[k];
1577                 a[k] = x;
1578                 _j1 = j + m;
1579                 k1 = k + m;
1580                 x = a[_j1];
1581                 a[_j1] = a[k1];
1582                 a[k1] = x;
1583                 _j1 += nh;
1584                 k1++;
1585                 x = a[_j1];
1586                 a[_j1] = a[k1];
1587                 a[k1] = x;
1588                 _j1 -= m;
1589                 k1 -= m;
1590                 x = a[_j1];
1591                 a[_j1] = a[k1];
1592                 a[k1] = x;
1593                 _j1++;
1594                 k1 += nh;
1595                 x = a[_j1];
1596                 a[_j1] = a[k1];
1597                 a[k1] = x;
1598                 _j1 += m;
1599                 k1 += m;
1600                 x = a[_j1];
1601                 a[_j1] = a[k1];
1602                 a[k1] = x;
1603                 _j1 -= nh;
1604                 k1--;
1605                 x = a[_j1];
1606                 a[_j1] = a[k1];
1607                 a[k1] = x;
1608                 _j1 -= m;
1609                 k1 -= m;
1610                 x = a[_j1];
1611                 a[_j1] = a[k1];
1612                 a[k1] = x;
1613                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
1614             }
1615             k1 = _j0 + k0;
1616             _j1 = k1 + 1;
1617             k1 += nh;
1618             x = a[_j1];
1619             a[_j1] = a[k1];
1620             a[k1] = x;
1621             _j1 += m;
1622             k1 += m;
1623             x = a[_j1];
1624             a[_j1] = a[k1];
1625             a[k1] = x;
1626         }
1627     } else {
1628         _j0 = 0;
1629         for (k0 = 2; k0 < m; k0 += 2) {
1630             for (i = nh >> 1; i > (_j0 ^= i); i >>= 1);
1631             k = k0;
1632             for (j = _j0; j < _j0 + k0; j += 2) {
1633                 x = a[j];
1634                 a[j] = a[k];
1635                 a[k] = x;
1636                 _j1 = j + nh;
1637                 k1 = k + 1;
1638                 x = a[_j1];
1639                 a[_j1] = a[k1];
1640                 a[k1] = x;
1641                 _j1++;
1642                 k1 += nh;
1643                 x = a[_j1];
1644                 a[_j1] = a[k1];
1645                 a[k1] = x;
1646                 _j1 -= nh;
1647                 k1--;
1648                 x = a[_j1];
1649                 a[_j1] = a[k1];
1650                 a[k1] = x;
1651                 for (i = nh >> 1; i > (k ^= i); i >>= 1);
1652             }
1653             k1 = _j0 + k0;
1654             _j1 = k1 + 1;
1655             k1 += nh;
1656             x = a[_j1];
1657             a[_j1] = a[k1];
1658             a[k1] = x;
1659         }
1660     }
1661 }
1662 
1663 
1664 void cftb1st(int n, double *a)
1665 {
1666     int i, i0, j, _j0, _j1, j2, j3, m, mh;
1667     double ew, w1r, w1i, wk1r, wk1i, wk3r, wk3i, 
1668         wd1r, wd1i, wd3r, wd3i, ss1, ss3;
1669     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1670     
1671     mh = n >> 3;
1672     m = 2 * mh;
1673     _j1 = m;
1674     j2 = _j1 + m;
1675     j3 = j2 + m;
1676     x0r = a[0] + a[j2];
1677     x0i = -a[1] - a[j2 + 1];
1678     x1r = a[0] - a[j2];
1679     x1i = -a[1] + a[j2 + 1];
1680     x2r = a[_j1] + a[j3];
1681     x2i = a[_j1 + 1] + a[j3 + 1];
1682     x3r = a[_j1] - a[j3];
1683     x3i = a[_j1 + 1] - a[j3 + 1];
1684     a[0] = x0r + x2r;
1685     a[1] = x0i - x2i;
1686     a[_j1] = x0r - x2r;
1687     a[_j1 + 1] = x0i + x2i;
1688     a[j2] = x1r + x3i;
1689     a[j2 + 1] = x1i + x3r;
1690     a[j3] = x1r - x3i;
1691     a[j3 + 1] = x1i - x3r;
1692     wd1r = 1;
1693     wd1i = 0;
1694     wd3r = 1;
1695     wd3i = 0;
1696     ew = M_PI_2 / m;
1697     w1r = cos(2 * ew);
1698     w1i = sin(2 * ew);
1699     wk1r = w1r;
1700     wk1i = w1i;
1701     ss1 = 2 * w1i;
1702     wk3i = 2 * ss1 * wk1r;
1703     wk3r = wk1r - wk3i * wk1i;
1704     wk3i = wk1i - wk3i * wk1r;
1705     ss3 = 2 * wk3i;
1706     i = 0;
1707     for (;;) {
1708         i0 = i + 4 * CDFT_LOOP_DIV;
1709         if (i0 > mh - 4) {
1710             i0 = mh - 4;
1711         }
1712         for (j = i + 2; j < i0; j += 4) {
1713             wd1r -= ss1 * wk1i;
1714             wd1i += ss1 * wk1r;
1715             wd3r -= ss3 * wk3i;
1716             wd3i += ss3 * wk3r;
1717             _j1 = j + m;
1718             j2 = _j1 + m;
1719             j3 = j2 + m;
1720             x0r = a[j] + a[j2];
1721             x0i = -a[j + 1] - a[j2 + 1];
1722             x1r = a[j] - a[j2];
1723             x1i = -a[j + 1] + a[j2 + 1];
1724             x2r = a[_j1] + a[j3];
1725             x2i = a[_j1 + 1] + a[j3 + 1];
1726             x3r = a[_j1] - a[j3];
1727             x3i = a[_j1 + 1] - a[j3 + 1];
1728             a[j] = x0r + x2r;
1729             a[j + 1] = x0i - x2i;
1730             a[_j1] = x0r - x2r;
1731             a[_j1 + 1] = x0i + x2i;
1732             x0r = x1r + x3i;
1733             x0i = x1i + x3r;
1734             a[j2] = wk1r * x0r - wk1i * x0i;
1735             a[j2 + 1] = wk1r * x0i + wk1i * x0r;
1736             x0r = x1r - x3i;
1737             x0i = x1i - x3r;
1738             a[j3] = wk3r * x0r + wk3i * x0i;
1739             a[j3 + 1] = wk3r * x0i - wk3i * x0r;
1740             x0r = a[j + 2] + a[j2 + 2];
1741             x0i = -a[j + 3] - a[j2 + 3];
1742             x1r = a[j + 2] - a[j2 + 2];
1743             x1i = -a[j + 3] + a[j2 + 3];
1744             x2r = a[_j1 + 2] + a[j3 + 2];
1745             x2i = a[_j1 + 3] + a[j3 + 3];
1746             x3r = a[_j1 + 2] - a[j3 + 2];
1747             x3i = a[_j1 + 3] - a[j3 + 3];
1748             a[j + 2] = x0r + x2r;
1749             a[j + 3] = x0i - x2i;
1750             a[_j1 + 2] = x0r - x2r;
1751             a[_j1 + 3] = x0i + x2i;
1752             x0r = x1r + x3i;
1753             x0i = x1i + x3r;
1754             a[j2 + 2] = wd1r * x0r - wd1i * x0i;
1755             a[j2 + 3] = wd1r * x0i + wd1i * x0r;
1756             x0r = x1r - x3i;
1757             x0i = x1i - x3r;
1758             a[j3 + 2] = wd3r * x0r + wd3i * x0i;
1759             a[j3 + 3] = wd3r * x0i - wd3i * x0r;
1760             _j0 = m - j;
1761             _j1 = _j0 + m;
1762             j2 = _j1 + m;
1763             j3 = j2 + m;
1764             x0r = a[_j0] + a[j2];
1765             x0i = -a[_j0 + 1] - a[j2 + 1];
1766             x1r = a[_j0] - a[j2];
1767             x1i = -a[_j0 + 1] + a[j2 + 1];
1768             x2r = a[_j1] + a[j3];
1769             x2i = a[_j1 + 1] + a[j3 + 1];
1770             x3r = a[_j1] - a[j3];
1771             x3i = a[_j1 + 1] - a[j3 + 1];
1772             a[_j0] = x0r + x2r;
1773             a[_j0 + 1] = x0i - x2i;
1774             a[_j1] = x0r - x2r;
1775             a[_j1 + 1] = x0i + x2i;
1776             x0r = x1r + x3i;
1777             x0i = x1i + x3r;
1778             a[j2] = wk1i * x0r - wk1r * x0i;
1779             a[j2 + 1] = wk1i * x0i + wk1r * x0r;
1780             x0r = x1r - x3i;
1781             x0i = x1i - x3r;
1782             a[j3] = wk3i * x0r + wk3r * x0i;
1783             a[j3 + 1] = wk3i * x0i - wk3r * x0r;
1784             x0r = a[_j0 - 2] + a[j2 - 2];
1785             x0i = -a[_j0 - 1] - a[j2 - 1];
1786             x1r = a[_j0 - 2] - a[j2 - 2];
1787             x1i = -a[_j0 - 1] + a[j2 - 1];
1788             x2r = a[_j1 - 2] + a[j3 - 2];
1789             x2i = a[_j1 - 1] + a[j3 - 1];
1790             x3r = a[_j1 - 2] - a[j3 - 2];
1791             x3i = a[_j1 - 1] - a[j3 - 1];
1792             a[_j0 - 2] = x0r + x2r;
1793             a[_j0 - 1] = x0i - x2i;
1794             a[_j1 - 2] = x0r - x2r;
1795             a[_j1 - 1] = x0i + x2i;
1796             x0r = x1r + x3i;
1797             x0i = x1i + x3r;
1798             a[j2 - 2] = wd1i * x0r - wd1r * x0i;
1799             a[j2 - 1] = wd1i * x0i + wd1r * x0r;
1800             x0r = x1r - x3i;
1801             x0i = x1i - x3r;
1802             a[j3 - 2] = wd3i * x0r + wd3r * x0i;
1803             a[j3 - 1] = wd3i * x0i - wd3r * x0r;
1804             wk1r -= ss1 * wd1i;
1805             wk1i += ss1 * wd1r;
1806             wk3r -= ss3 * wd3i;
1807             wk3i += ss3 * wd3r;
1808         }
1809         if (i0 == mh - 4) {
1810             break;
1811         }
1812         wd1r = cos(ew * i0);
1813         wd1i = sin(ew * i0);
1814         wd3i = 4 * wd1i * wd1r;
1815         wd3r = wd1r - wd3i * wd1i;
1816         wd3i = wd1i - wd3i * wd1r;
1817         wk1r = w1r * wd1r - w1i * wd1i;
1818         wk1i = w1r * wd1i + w1i * wd1r;
1819         wk3i = 4 * wk1i * wk1r;
1820         wk3r = wk1r - wk3i * wk1i;
1821         wk3i = wk1i - wk3i * wk1r;
1822         i = i0;
1823     }
1824     wd1r = WR5000;
1825     _j0 = mh;
1826     _j1 = _j0 + m;
1827     j2 = _j1 + m;
1828     j3 = j2 + m;
1829     x0r = a[_j0 - 2] + a[j2 - 2];
1830     x0i = -a[_j0 - 1] - a[j2 - 1];
1831     x1r = a[_j0 - 2] - a[j2 - 2];
1832     x1i = -a[_j0 - 1] + a[j2 - 1];
1833     x2r = a[_j1 - 2] + a[j3 - 2];
1834     x2i = a[_j1 - 1] + a[j3 - 1];
1835     x3r = a[_j1 - 2] - a[j3 - 2];
1836     x3i = a[_j1 - 1] - a[j3 - 1];
1837     a[_j0 - 2] = x0r + x2r;
1838     a[_j0 - 1] = x0i - x2i;
1839     a[_j1 - 2] = x0r - x2r;
1840     a[_j1 - 1] = x0i + x2i;
1841     x0r = x1r + x3i;
1842     x0i = x1i + x3r;
1843     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
1844     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
1845     x0r = x1r - x3i;
1846     x0i = x1i - x3r;
1847     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
1848     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
1849     x0r = a[_j0] + a[j2];
1850     x0i = -a[_j0 + 1] - a[j2 + 1];
1851     x1r = a[_j0] - a[j2];
1852     x1i = -a[_j0 + 1] + a[j2 + 1];
1853     x2r = a[_j1] + a[j3];
1854     x2i = a[_j1 + 1] + a[j3 + 1];
1855     x3r = a[_j1] - a[j3];
1856     x3i = a[_j1 + 1] - a[j3 + 1];
1857     a[_j0] = x0r + x2r;
1858     a[_j0 + 1] = x0i - x2i;
1859     a[_j1] = x0r - x2r;
1860     a[_j1 + 1] = x0i + x2i;
1861     x0r = x1r + x3i;
1862     x0i = x1i + x3r;
1863     a[j2] = wd1r * (x0r - x0i);
1864     a[j2 + 1] = wd1r * (x0i + x0r);
1865     x0r = x1r - x3i;
1866     x0i = x1i - x3r;
1867     a[j3] = -wd1r * (x0r + x0i);
1868     a[j3 + 1] = -wd1r * (x0i - x0r);
1869     x0r = a[_j0 + 2] + a[j2 + 2];
1870     x0i = -a[_j0 + 3] - a[j2 + 3];
1871     x1r = a[_j0 + 2] - a[j2 + 2];
1872     x1i = -a[_j0 + 3] + a[j2 + 3];
1873     x2r = a[_j1 + 2] + a[j3 + 2];
1874     x2i = a[_j1 + 3] + a[j3 + 3];
1875     x3r = a[_j1 + 2] - a[j3 + 2];
1876     x3i = a[_j1 + 3] - a[j3 + 3];
1877     a[_j0 + 2] = x0r + x2r;
1878     a[_j0 + 3] = x0i - x2i;
1879     a[_j1 + 2] = x0r - x2r;
1880     a[_j1 + 3] = x0i + x2i;
1881     x0r = x1r + x3i;
1882     x0i = x1i + x3r;
1883     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
1884     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
1885     x0r = x1r - x3i;
1886     x0i = x1i - x3r;
1887     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
1888     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
1889 }
1890 
1891 
1892 #ifdef USE_CDFT_THREADS
1893 struct cdft_arg_st {
1894     int n0;
1895     int n;
1896     double *a;
1897 };
1898 typedef struct cdft_arg_st cdft_arg_t;
1899 
1900 
1901 void cftrec4_th(int n, double *a)
1902 {
1903     void *cftrec1_th(void *p);
1904     void *cftrec2_th(void *p);
1905     int i, idiv4, m, nthread;
1906     cdft_thread_t th[4];
1907     cdft_arg_t ag[4];
1908     
1909     nthread = 2;
1910     idiv4 = 0;
1911     m = n >> 1;
1912     if (n > CDFT_4THREADS_BEGIN_N) {
1913         nthread = 4;
1914         idiv4 = 1;
1915         m >>= 1;
1916     }
1917     for (i = 0; i < nthread; i++) {
1918         ag[i].n0 = n;
1919         ag[i].n = m;
1920         ag[i].a = &a[i * m];
1921         if (i != idiv4) {
1922             cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
1923         } else {
1924             cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
1925         }
1926     }
1927     for (i = 0; i < nthread; i++) {
1928         cdft_thread_wait(th[i]);
1929     }
1930 }
1931 
1932 
1933 void *cftrec1_th(void *p)
1934 {
1935     int cfttree(int n, int j, int k, double *a);
1936     void cftleaf(int n, int isplt, double *a);
1937     void cftmdl1(int n, double *a);
1938     int isplt, j, k, m, n, n0;
1939     double *a;
1940     
1941     n0 = ((cdft_arg_t *) p)->n0;
1942     n = ((cdft_arg_t *) p)->n;
1943     a = ((cdft_arg_t *) p)->a;
1944     m = n0;
1945     while (m > 512) {
1946         m >>= 2;
1947         cftmdl1(m, &a[n - m]);
1948     }
1949     cftleaf(m, 1, &a[n - m]);
1950     k = 0;
1951     for (j = n - m; j > 0; j -= m) {
1952         k++;
1953         isplt = cfttree(m, j, k, a);
1954         cftleaf(m, isplt, &a[j - m]);
1955     }
1956     return (void *) 0;
1957 }
1958 
1959 
1960 void *cftrec2_th(void *p)
1961 {
1962     int cfttree(int n, int j, int k, double *a);
1963     void cftleaf(int n, int isplt, double *a);
1964     void cftmdl2(int n, double *a);
1965     int isplt, j, k, m, n, n0;
1966     double *a;
1967     
1968     n0 = ((cdft_arg_t *) p)->n0;
1969     n = ((cdft_arg_t *) p)->n;
1970     a = ((cdft_arg_t *) p)->a;
1971     k = 1;
1972     m = n0;
1973     while (m > 512) {
1974         m >>= 2;
1975         k <<= 2;
1976         cftmdl2(m, &a[n - m]);
1977     }
1978     cftleaf(m, 0, &a[n - m]);
1979     k >>= 1;
1980     for (j = n - m; j > 0; j -= m) {
1981         k++;
1982         isplt = cfttree(m, j, k, a);
1983         cftleaf(m, isplt, &a[j - m]);
1984     }
1985     return (void *) 0;
1986 }
1987 #endif /* USE_CDFT_THREADS */
1988 
1989 
1990 void cftrec4(int n, double *a)
1991 {
1992     int isplt, j, k, m;
1993     
1994     m = n;
1995     while (m > 512) {
1996         m >>= 2;
1997         cftmdl1(m, &a[n - m]);
1998     }
1999     cftleaf(m, 1, &a[n - m]);
2000     k = 0;
2001     for (j = n - m; j > 0; j -= m) {
2002         k++;
2003         isplt = cfttree(m, j, k, a);
2004         cftleaf(m, isplt, &a[j - m]);
2005     }
2006 }
2007 
2008 
2009 int cfttree(int n, int j, int k, double *a)
2010 {
2011     int i, isplt, m;
2012     
2013     if ((k & 3) != 0) {
2014         isplt = k & 1;
2015         if (isplt != 0) {
2016             cftmdl1(n, &a[j - n]);
2017         } else {
2018             cftmdl2(n, &a[j - n]);
2019         }
2020     } else {
2021         m = n;
2022         for (i = k; (i & 3) == 0; i >>= 2) {
2023             m <<= 2;
2024         }
2025         isplt = i & 1;
2026         if (isplt != 0) {
2027             while (m > 128) {
2028                 cftmdl1(m, &a[j - m]);
2029                 m >>= 2;
2030             }
2031         } else {
2032             while (m > 128) {
2033                 cftmdl2(m, &a[j - m]);
2034                 m >>= 2;
2035             }
2036         }
2037     }
2038     return isplt;
2039 }
2040 
2041 
2042 void cftleaf(int n, int isplt, double *a)
2043 {
2044     if (n == 512) {
2045         cftmdl1(128, a);
2046         cftf161(a);
2047         cftf162(&a[32]);
2048         cftf161(&a[64]);
2049         cftf161(&a[96]);
2050         cftmdl2(128, &a[128]);
2051         cftf161(&a[128]);
2052         cftf162(&a[160]);
2053         cftf161(&a[192]);
2054         cftf162(&a[224]);
2055         cftmdl1(128, &a[256]);
2056         cftf161(&a[256]);
2057         cftf162(&a[288]);
2058         cftf161(&a[320]);
2059         cftf161(&a[352]);
2060         if (isplt != 0) {
2061             cftmdl1(128, &a[384]);
2062             cftf161(&a[480]);
2063         } else {
2064             cftmdl2(128, &a[384]);
2065             cftf162(&a[480]);
2066         }
2067         cftf161(&a[384]);
2068         cftf162(&a[416]);
2069         cftf161(&a[448]);
2070     } else {
2071         cftmdl1(64, a);
2072         cftf081(a);
2073         cftf082(&a[16]);
2074         cftf081(&a[32]);
2075         cftf081(&a[48]);
2076         cftmdl2(64, &a[64]);
2077         cftf081(&a[64]);
2078         cftf082(&a[80]);
2079         cftf081(&a[96]);
2080         cftf082(&a[112]);
2081         cftmdl1(64, &a[128]);
2082         cftf081(&a[128]);
2083         cftf082(&a[144]);
2084         cftf081(&a[160]);
2085         cftf081(&a[176]);
2086         if (isplt != 0) {
2087             cftmdl1(64, &a[192]);
2088             cftf081(&a[240]);
2089         } else {
2090             cftmdl2(64, &a[192]);
2091             cftf082(&a[240]);
2092         }
2093         cftf081(&a[192]);
2094         cftf082(&a[208]);
2095         cftf081(&a[224]);
2096     }
2097 }
2098 
2099 
2100 void cftmdl1(int n, double *a)
2101 {
2102     int i, i0, j, _j0, _j1, j2, j3, m, mh;
2103     double ew, w1r, w1i, wk1r, wk1i, wk3r, wk3i, 
2104         wd1r, wd1i, wd3r, wd3i, ss1, ss3;
2105     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2106     
2107     mh = n >> 3;
2108     m = 2 * mh;
2109     _j1 = m;
2110     j2 = _j1 + m;
2111     j3 = j2 + m;
2112     x0r = a[0] + a[j2];
2113     x0i = a[1] + a[j2 + 1];
2114     x1r = a[0] - a[j2];
2115     x1i = a[1] - a[j2 + 1];
2116     x2r = a[_j1] + a[j3];
2117     x2i = a[_j1 + 1] + a[j3 + 1];
2118     x3r = a[_j1] - a[j3];
2119     x3i = a[_j1 + 1] - a[j3 + 1];
2120     a[0] = x0r + x2r;
2121     a[1] = x0i + x2i;
2122     a[_j1] = x0r - x2r;
2123     a[_j1 + 1] = x0i - x2i;
2124     a[j2] = x1r - x3i;
2125     a[j2 + 1] = x1i + x3r;
2126     a[j3] = x1r + x3i;
2127     a[j3 + 1] = x1i - x3r;
2128     wd1r = 1;
2129     wd1i = 0;
2130     wd3r = 1;
2131     wd3i = 0;
2132     ew = M_PI_2 / m;
2133     w1r = cos(2 * ew);
2134     w1i = sin(2 * ew);
2135     wk1r = w1r;
2136     wk1i = w1i;
2137     ss1 = 2 * w1i;
2138     wk3i = 2 * ss1 * wk1r;
2139     wk3r = wk1r - wk3i * wk1i;
2140     wk3i = wk1i - wk3i * wk1r;
2141     ss3 = 2 * wk3i;
2142     i = 0;
2143     for (;;) {
2144         i0 = i + 4 * CDFT_LOOP_DIV;
2145         if (i0 > mh - 4) {
2146             i0 = mh - 4;
2147         }
2148         for (j = i + 2; j < i0; j += 4) {
2149             wd1r -= ss1 * wk1i;
2150             wd1i += ss1 * wk1r;
2151             wd3r -= ss3 * wk3i;
2152             wd3i += ss3 * wk3r;
2153             _j1 = j + m;
2154             j2 = _j1 + m;
2155             j3 = j2 + m;
2156             x0r = a[j] + a[j2];
2157             x0i = a[j + 1] + a[j2 + 1];
2158             x1r = a[j] - a[j2];
2159             x1i = a[j + 1] - a[j2 + 1];
2160             x2r = a[_j1] + a[j3];
2161             x2i = a[_j1 + 1] + a[j3 + 1];
2162             x3r = a[_j1] - a[j3];
2163             x3i = a[_j1 + 1] - a[j3 + 1];
2164             a[j] = x0r + x2r;
2165             a[j + 1] = x0i + x2i;
2166             a[_j1] = x0r - x2r;
2167             a[_j1 + 1] = x0i - x2i;
2168             x0r = x1r - x3i;
2169             x0i = x1i + x3r;
2170             a[j2] = wk1r * x0r - wk1i * x0i;
2171             a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2172             x0r = x1r + x3i;
2173             x0i = x1i - x3r;
2174             a[j3] = wk3r * x0r + wk3i * x0i;
2175             a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2176             x0r = a[j + 2] + a[j2 + 2];
2177             x0i = a[j + 3] + a[j2 + 3];
2178             x1r = a[j + 2] - a[j2 + 2];
2179             x1i = a[j + 3] - a[j2 + 3];
2180             x2r = a[_j1 + 2] + a[j3 + 2];
2181             x2i = a[_j1 + 3] + a[j3 + 3];
2182             x3r = a[_j1 + 2] - a[j3 + 2];
2183             x3i = a[_j1 + 3] - a[j3 + 3];
2184             a[j + 2] = x0r + x2r;
2185             a[j + 3] = x0i + x2i;
2186             a[_j1 + 2] = x0r - x2r;
2187             a[_j1 + 3] = x0i - x2i;
2188             x0r = x1r - x3i;
2189             x0i = x1i + x3r;
2190             a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2191             a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2192             x0r = x1r + x3i;
2193             x0i = x1i - x3r;
2194             a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2195             a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2196             _j0 = m - j;
2197             _j1 = _j0 + m;
2198             j2 = _j1 + m;
2199             j3 = j2 + m;
2200             x0r = a[_j0] + a[j2];
2201             x0i = a[_j0 + 1] + a[j2 + 1];
2202             x1r = a[_j0] - a[j2];
2203             x1i = a[_j0 + 1] - a[j2 + 1];
2204             x2r = a[_j1] + a[j3];
2205             x2i = a[_j1 + 1] + a[j3 + 1];
2206             x3r = a[_j1] - a[j3];
2207             x3i = a[_j1 + 1] - a[j3 + 1];
2208             a[_j0] = x0r + x2r;
2209             a[_j0 + 1] = x0i + x2i;
2210             a[_j1] = x0r - x2r;
2211             a[_j1 + 1] = x0i - x2i;
2212             x0r = x1r - x3i;
2213             x0i = x1i + x3r;
2214             a[j2] = wk1i * x0r - wk1r * x0i;
2215             a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2216             x0r = x1r + x3i;
2217             x0i = x1i - x3r;
2218             a[j3] = wk3i * x0r + wk3r * x0i;
2219             a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2220             x0r = a[_j0 - 2] + a[j2 - 2];
2221             x0i = a[_j0 - 1] + a[j2 - 1];
2222             x1r = a[_j0 - 2] - a[j2 - 2];
2223             x1i = a[_j0 - 1] - a[j2 - 1];
2224             x2r = a[_j1 - 2] + a[j3 - 2];
2225             x2i = a[_j1 - 1] + a[j3 - 1];
2226             x3r = a[_j1 - 2] - a[j3 - 2];
2227             x3i = a[_j1 - 1] - a[j3 - 1];
2228             a[_j0 - 2] = x0r + x2r;
2229             a[_j0 - 1] = x0i + x2i;
2230             a[_j1 - 2] = x0r - x2r;
2231             a[_j1 - 1] = x0i - x2i;
2232             x0r = x1r - x3i;
2233             x0i = x1i + x3r;
2234             a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2235             a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2236             x0r = x1r + x3i;
2237             x0i = x1i - x3r;
2238             a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2239             a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2240             wk1r -= ss1 * wd1i;
2241             wk1i += ss1 * wd1r;
2242             wk3r -= ss3 * wd3i;
2243             wk3i += ss3 * wd3r;
2244         }
2245         if (i0 == mh - 4) {
2246             break;
2247         }
2248         wd1r = cos(ew * i0);
2249         wd1i = sin(ew * i0);
2250         wd3i = 4 * wd1i * wd1r;
2251         wd3r = wd1r - wd3i * wd1i;
2252         wd3i = wd1i - wd3i * wd1r;
2253         wk1r = w1r * wd1r - w1i * wd1i;
2254         wk1i = w1r * wd1i + w1i * wd1r;
2255         wk3i = 4 * wk1i * wk1r;
2256         wk3r = wk1r - wk3i * wk1i;
2257         wk3i = wk1i - wk3i * wk1r;
2258         i = i0;
2259     }
2260     wd1r = WR5000;
2261     _j0 = mh;
2262     _j1 = _j0 + m;
2263     j2 = _j1 + m;
2264     j3 = j2 + m;
2265     x0r = a[_j0 - 2] + a[j2 - 2];
2266     x0i = a[_j0 - 1] + a[j2 - 1];
2267     x1r = a[_j0 - 2] - a[j2 - 2];
2268     x1i = a[_j0 - 1] - a[j2 - 1];
2269     x2r = a[_j1 - 2] + a[j3 - 2];
2270     x2i = a[_j1 - 1] + a[j3 - 1];
2271     x3r = a[_j1 - 2] - a[j3 - 2];
2272     x3i = a[_j1 - 1] - a[j3 - 1];
2273     a[_j0 - 2] = x0r + x2r;
2274     a[_j0 - 1] = x0i + x2i;
2275     a[_j1 - 2] = x0r - x2r;
2276     a[_j1 - 1] = x0i - x2i;
2277     x0r = x1r - x3i;
2278     x0i = x1i + x3r;
2279     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2280     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2281     x0r = x1r + x3i;
2282     x0i = x1i - x3r;
2283     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2284     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2285     x0r = a[_j0] + a[j2];
2286     x0i = a[_j0 + 1] + a[j2 + 1];
2287     x1r = a[_j0] - a[j2];
2288     x1i = a[_j0 + 1] - a[j2 + 1];
2289     x2r = a[_j1] + a[j3];
2290     x2i = a[_j1 + 1] + a[j3 + 1];
2291     x3r = a[_j1] - a[j3];
2292     x3i = a[_j1 + 1] - a[j3 + 1];
2293     a[_j0] = x0r + x2r;
2294     a[_j0 + 1] = x0i + x2i;
2295     a[_j1] = x0r - x2r;
2296     a[_j1 + 1] = x0i - x2i;
2297     x0r = x1r - x3i;
2298     x0i = x1i + x3r;
2299     a[j2] = wd1r * (x0r - x0i);
2300     a[j2 + 1] = wd1r * (x0i + x0r);
2301     x0r = x1r + x3i;
2302     x0i = x1i - x3r;
2303     a[j3] = -wd1r * (x0r + x0i);
2304     a[j3 + 1] = -wd1r * (x0i - x0r);
2305     x0r = a[_j0 + 2] + a[j2 + 2];
2306     x0i = a[_j0 + 3] + a[j2 + 3];
2307     x1r = a[_j0 + 2] - a[j2 + 2];
2308     x1i = a[_j0 + 3] - a[j2 + 3];
2309     x2r = a[_j1 + 2] + a[j3 + 2];
2310     x2i = a[_j1 + 3] + a[j3 + 3];
2311     x3r = a[_j1 + 2] - a[j3 + 2];
2312     x3i = a[_j1 + 3] - a[j3 + 3];
2313     a[_j0 + 2] = x0r + x2r;
2314     a[_j0 + 3] = x0i + x2i;
2315     a[_j1 + 2] = x0r - x2r;
2316     a[_j1 + 3] = x0i - x2i;
2317     x0r = x1r - x3i;
2318     x0i = x1i + x3r;
2319     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2320     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2321     x0r = x1r + x3i;
2322     x0i = x1i - x3r;
2323     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2324     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2325 }
2326 
2327 
2328 void cftmdl2(int n, double *a)
2329 {
2330     int i, i0, j, _j0, _j1, j2, j3, m, mh;
2331     double ew, w1r, w1i, wn4r, wk1r, wk1i, wk3r, wk3i, 
2332         wl1r, wl1i, wl3r, wl3i, wd1r, wd1i, wd3r, wd3i, 
2333         we1r, we1i, we3r, we3i, ss1, ss3;
2334     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2335     
2336     mh = n >> 3;
2337     m = 2 * mh;
2338     wn4r = WR5000;
2339     _j1 = m;
2340     j2 = _j1 + m;
2341     j3 = j2 + m;
2342     x0r = a[0] - a[j2 + 1];
2343     x0i = a[1] + a[j2];
2344     x1r = a[0] + a[j2 + 1];
2345     x1i = a[1] - a[j2];
2346     x2r = a[_j1] - a[j3 + 1];
2347     x2i = a[_j1 + 1] + a[j3];
2348     x3r = a[_j1] + a[j3 + 1];
2349     x3i = a[_j1 + 1] - a[j3];
2350     y0r = wn4r * (x2r - x2i);
2351     y0i = wn4r * (x2i + x2r);
2352     a[0] = x0r + y0r;
2353     a[1] = x0i + y0i;
2354     a[_j1] = x0r - y0r;
2355     a[_j1 + 1] = x0i - y0i;
2356     y0r = wn4r * (x3r - x3i);
2357     y0i = wn4r * (x3i + x3r);
2358     a[j2] = x1r - y0i;
2359     a[j2 + 1] = x1i + y0r;
2360     a[j3] = x1r + y0i;
2361     a[j3 + 1] = x1i - y0r;
2362     wl1r = 1;
2363     wl1i = 0;
2364     wl3r = 1;
2365     wl3i = 0;
2366     we1r = wn4r;
2367     we1i = wn4r;
2368     we3r = -wn4r;
2369     we3i = -wn4r;
2370     ew = M_PI_2 / (2 * m);
2371     w1r = cos(2 * ew);
2372     w1i = sin(2 * ew);
2373     wk1r = w1r;
2374     wk1i = w1i;
2375     wd1r = wn4r * (w1r - w1i);
2376     wd1i = wn4r * (w1i + w1r);
2377     ss1 = 2 * w1i;
2378     wk3i = 2 * ss1 * wk1r;
2379     wk3r = wk1r - wk3i * wk1i;
2380     wk3i = wk1i - wk3i * wk1r;
2381     ss3 = 2 * wk3i;
2382     wd3r = -wn4r * (wk3r - wk3i);
2383     wd3i = -wn4r * (wk3i + wk3r);
2384     i = 0;
2385     for (;;) {
2386         i0 = i + 4 * CDFT_LOOP_DIV;
2387         if (i0 > mh - 4) {
2388             i0 = mh - 4;
2389         }
2390         for (j = i + 2; j < i0; j += 4) {
2391             wl1r -= ss1 * wk1i;
2392             wl1i += ss1 * wk1r;
2393             wl3r -= ss3 * wk3i;
2394             wl3i += ss3 * wk3r;
2395             we1r -= ss1 * wd1i;
2396             we1i += ss1 * wd1r;
2397             we3r -= ss3 * wd3i;
2398             we3i += ss3 * wd3r;
2399             _j1 = j + m;
2400             j2 = _j1 + m;
2401             j3 = j2 + m;
2402             x0r = a[j] - a[j2 + 1];
2403             x0i = a[j + 1] + a[j2];
2404             x1r = a[j] + a[j2 + 1];
2405             x1i = a[j + 1] - a[j2];
2406             x2r = a[_j1] - a[j3 + 1];
2407             x2i = a[_j1 + 1] + a[j3];
2408             x3r = a[_j1] + a[j3 + 1];
2409             x3i = a[_j1 + 1] - a[j3];
2410             y0r = wk1r * x0r - wk1i * x0i;
2411             y0i = wk1r * x0i + wk1i * x0r;
2412             y2r = wd1r * x2r - wd1i * x2i;
2413             y2i = wd1r * x2i + wd1i * x2r;
2414             a[j] = y0r + y2r;
2415             a[j + 1] = y0i + y2i;
2416             a[_j1] = y0r - y2r;
2417             a[_j1 + 1] = y0i - y2i;
2418             y0r = wk3r * x1r + wk3i * x1i;
2419             y0i = wk3r * x1i - wk3i * x1r;
2420             y2r = wd3r * x3r + wd3i * x3i;
2421             y2i = wd3r * x3i - wd3i * x3r;
2422             a[j2] = y0r + y2r;
2423             a[j2 + 1] = y0i + y2i;
2424             a[j3] = y0r - y2r;
2425             a[j3 + 1] = y0i - y2i;
2426             x0r = a[j + 2] - a[j2 + 3];
2427             x0i = a[j + 3] + a[j2 + 2];
2428             x1r = a[j + 2] + a[j2 + 3];
2429             x1i = a[j + 3] - a[j2 + 2];
2430             x2r = a[_j1 + 2] - a[j3 + 3];
2431             x2i = a[_j1 + 3] + a[j3 + 2];
2432             x3r = a[_j1 + 2] + a[j3 + 3];
2433             x3i = a[_j1 + 3] - a[j3 + 2];
2434             y0r = wl1r * x0r - wl1i * x0i;
2435             y0i = wl1r * x0i + wl1i * x0r;
2436             y2r = we1r * x2r - we1i * x2i;
2437             y2i = we1r * x2i + we1i * x2r;
2438             a[j + 2] = y0r + y2r;
2439             a[j + 3] = y0i + y2i;
2440             a[_j1 + 2] = y0r - y2r;
2441             a[_j1 + 3] = y0i - y2i;
2442             y0r = wl3r * x1r + wl3i * x1i;
2443             y0i = wl3r * x1i - wl3i * x1r;
2444             y2r = we3r * x3r + we3i * x3i;
2445             y2i = we3r * x3i - we3i * x3r;
2446             a[j2 + 2] = y0r + y2r;
2447             a[j2 + 3] = y0i + y2i;
2448             a[j3 + 2] = y0r - y2r;
2449             a[j3 + 3] = y0i - y2i;
2450             _j0 = m - j;
2451             _j1 = _j0 + m;
2452             j2 = _j1 + m;
2453             j3 = j2 + m;
2454             x0r = a[_j0] - a[j2 + 1];
2455             x0i = a[_j0 + 1] + a[j2];
2456             x1r = a[_j0] + a[j2 + 1];
2457             x1i = a[_j0 + 1] - a[j2];
2458             x2r = a[_j1] - a[j3 + 1];
2459             x2i = a[_j1 + 1] + a[j3];
2460             x3r = a[_j1] + a[j3 + 1];
2461             x3i = a[_j1 + 1] - a[j3];
2462             y0r = wd1i * x0r - wd1r * x0i;
2463             y0i = wd1i * x0i + wd1r * x0r;
2464             y2r = wk1i * x2r - wk1r * x2i;
2465             y2i = wk1i * x2i + wk1r * x2r;
2466             a[_j0] = y0r + y2r;
2467             a[_j0 + 1] = y0i + y2i;
2468             a[_j1] = y0r - y2r;
2469             a[_j1 + 1] = y0i - y2i;
2470             y0r = wd3i * x1r + wd3r * x1i;
2471             y0i = wd3i * x1i - wd3r * x1r;
2472             y2r = wk3i * x3r + wk3r * x3i;
2473             y2i = wk3i * x3i - wk3r * x3r;
2474             a[j2] = y0r + y2r;
2475             a[j2 + 1] = y0i + y2i;
2476             a[j3] = y0r - y2r;
2477             a[j3 + 1] = y0i - y2i;
2478             x0r = a[_j0 - 2] - a[j2 - 1];
2479             x0i = a[_j0 - 1] + a[j2 - 2];
2480             x1r = a[_j0 - 2] + a[j2 - 1];
2481             x1i = a[_j0 - 1] - a[j2 - 2];
2482             x2r = a[_j1 - 2] - a[j3 - 1];
2483             x2i = a[_j1 - 1] + a[j3 - 2];
2484             x3r = a[_j1 - 2] + a[j3 - 1];
2485             x3i = a[_j1 - 1] - a[j3 - 2];
2486             y0r = we1i * x0r - we1r * x0i;
2487             y0i = we1i * x0i + we1r * x0r;
2488             y2r = wl1i * x2r - wl1r * x2i;
2489             y2i = wl1i * x2i + wl1r * x2r;
2490             a[_j0 - 2] = y0r + y2r;
2491             a[_j0 - 1] = y0i + y2i;
2492             a[_j1 - 2] = y0r - y2r;
2493             a[_j1 - 1] = y0i - y2i;
2494             y0r = we3i * x1r + we3r * x1i;
2495             y0i = we3i * x1i - we3r * x1r;
2496             y2r = wl3i * x3r + wl3r * x3i;
2497             y2i = wl3i * x3i - wl3r * x3r;
2498             a[j2 - 2] = y0r + y2r;
2499             a[j2 - 1] = y0i + y2i;
2500             a[j3 - 2] = y0r - y2r;
2501             a[j3 - 1] = y0i - y2i;
2502             wk1r -= ss1 * wl1i;
2503             wk1i += ss1 * wl1r;
2504             wk3r -= ss3 * wl3i;
2505             wk3i += ss3 * wl3r;
2506             wd1r -= ss1 * we1i;
2507             wd1i += ss1 * we1r;
2508             wd3r -= ss3 * we3i;
2509             wd3i += ss3 * we3r;
2510         }
2511         if (i0 == mh - 4) {
2512             break;
2513         }
2514         wl1r = cos(ew * i0);
2515         wl1i = sin(ew * i0);
2516         wl3i = 4 * wl1i * wl1r;
2517         wl3r = wl1r - wl3i * wl1i;
2518         wl3i = wl1i - wl3i * wl1r;
2519         we1r = wn4r * (wl1r - wl1i);
2520         we1i = wn4r * (wl1i + wl1r);
2521         we3r = -wn4r * (wl3r - wl3i);
2522         we3i = -wn4r * (wl3i + wl3r);
2523         wk1r = w1r * wl1r - w1i * wl1i;
2524         wk1i = w1r * wl1i + w1i * wl1r;
2525         wk3i = 4 * wk1i * wk1r;
2526         wk3r = wk1r - wk3i * wk1i;
2527         wk3i = wk1i - wk3i * wk1r;
2528         wd1r = wn4r * (wk1r - wk1i);
2529         wd1i = wn4r * (wk1i + wk1r);
2530         wd3r = -wn4r * (wk3r - wk3i);
2531         wd3i = -wn4r * (wk3i + wk3r);
2532         i = i0;
2533     }
2534     wl1r = WR2500;
2535     wl1i = WI2500;
2536     _j0 = mh;
2537     _j1 = _j0 + m;
2538     j2 = _j1 + m;
2539     j3 = j2 + m;
2540     x0r = a[_j0 - 2] - a[j2 - 1];
2541     x0i = a[_j0 - 1] + a[j2 - 2];
2542     x1r = a[_j0 - 2] + a[j2 - 1];
2543     x1i = a[_j0 - 1] - a[j2 - 2];
2544     x2r = a[_j1 - 2] - a[j3 - 1];
2545     x2i = a[_j1 - 1] + a[j3 - 2];
2546     x3r = a[_j1 - 2] + a[j3 - 1];
2547     x3i = a[_j1 - 1] - a[j3 - 2];
2548     y0r = wk1r * x0r - wk1i * x0i;
2549     y0i = wk1r * x0i + wk1i * x0r;
2550     y2r = wd1r * x2r - wd1i * x2i;
2551     y2i = wd1r * x2i + wd1i * x2r;
2552     a[_j0 - 2] = y0r + y2r;
2553     a[_j0 - 1] = y0i + y2i;
2554     a[_j1 - 2] = y0r - y2r;
2555     a[_j1 - 1] = y0i - y2i;
2556     y0r = wk3r * x1r + wk3i * x1i;
2557     y0i = wk3r * x1i - wk3i * x1r;
2558     y2r = wd3r * x3r + wd3i * x3i;
2559     y2i = wd3r * x3i - wd3i * x3r;
2560     a[j2 - 2] = y0r + y2r;
2561     a[j2 - 1] = y0i + y2i;
2562     a[j3 - 2] = y0r - y2r;
2563     a[j3 - 1] = y0i - y2i;
2564     x0r = a[_j0] - a[j2 + 1];
2565     x0i = a[_j0 + 1] + a[j2];
2566     x1r = a[_j0] + a[j2 + 1];
2567     x1i = a[_j0 + 1] - a[j2];
2568     x2r = a[_j1] - a[j3 + 1];
2569     x2i = a[_j1 + 1] + a[j3];
2570     x3r = a[_j1] + a[j3 + 1];
2571     x3i = a[_j1 + 1] - a[j3];
2572     y0r = wl1r * x0r - wl1i * x0i;
2573     y0i = wl1r * x0i + wl1i * x0r;
2574     y2r = wl1i * x2r - wl1r * x2i;
2575     y2i = wl1i * x2i + wl1r * x2r;
2576     a[_j0] = y0r + y2r;
2577     a[_j0 + 1] = y0i + y2i;
2578     a[_j1] = y0r - y2r;
2579     a[_j1 + 1] = y0i - y2i;
2580     y0r = wl1i * x1r - wl1r * x1i;
2581     y0i = wl1i * x1i + wl1r * x1r;
2582     y2r = wl1r * x3r - wl1i * x3i;
2583     y2i = wl1r * x3i + wl1i * x3r;
2584     a[j2] = y0r - y2r;
2585     a[j2 + 1] = y0i - y2i;
2586     a[j3] = y0r + y2r;
2587     a[j3 + 1] = y0i + y2i;
2588     x0r = a[_j0 + 2] - a[j2 + 3];
2589     x0i = a[_j0 + 3] + a[j2 + 2];
2590     x1r = a[_j0 + 2] + a[j2 + 3];
2591     x1i = a[_j0 + 3] - a[j2 + 2];
2592     x2r = a[_j1 + 2] - a[j3 + 3];
2593     x2i = a[_j1 + 3] + a[j3 + 2];
2594     x3r = a[_j1 + 2] + a[j3 + 3];
2595     x3i = a[_j1 + 3] - a[j3 + 2];
2596     y0r = wd1i * x0r - wd1r * x0i;
2597     y0i = wd1i * x0i + wd1r * x0r;
2598     y2r = wk1i * x2r - wk1r * x2i;
2599     y2i = wk1i * x2i + wk1r * x2r;
2600     a[_j0 + 2] = y0r + y2r;
2601     a[_j0 + 3] = y0i + y2i;
2602     a[_j1 + 2] = y0r - y2r;
2603     a[_j1 + 3] = y0i - y2i;
2604     y0r = wd3i * x1r + wd3r * x1i;
2605     y0i = wd3i * x1i - wd3r * x1r;
2606     y2r = wk3i * x3r + wk3r * x3i;
2607     y2i = wk3i * x3i - wk3r * x3r;
2608     a[j2 + 2] = y0r + y2r;
2609     a[j2 + 3] = y0i + y2i;
2610     a[j3 + 2] = y0r - y2r;
2611     a[j3 + 3] = y0i - y2i;
2612 }
2613 
2614 
2615 void cftfx41(int n, double *a)
2616 {
2617     if (n == 128) {
2618         cftf161(a);
2619         cftf162(&a[32]);
2620         cftf161(&a[64]);
2621         cftf161(&a[96]);
2622     } else {
2623         cftf081(a);
2624         cftf082(&a[16]);
2625         cftf081(&a[32]);
2626         cftf081(&a[48]);
2627     }
2628 }
2629 
2630 
2631 void cftf161(double *a)
2632 {
2633     double wn4r, wk1r, wk1i, 
2634         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
2635         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
2636         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
2637         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
2638         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2639     
2640     wn4r = WR5000;
2641     wk1r = WR2500;
2642     wk1i = WI2500;
2643     x0r = a[0] + a[16];
2644     x0i = a[1] + a[17];
2645     x1r = a[0] - a[16];
2646     x1i = a[1] - a[17];
2647     x2r = a[8] + a[24];
2648     x2i = a[9] + a[25];
2649     x3r = a[8] - a[24];
2650     x3i = a[9] - a[25];
2651     y0r = x0r + x2r;
2652     y0i = x0i + x2i;
2653     y4r = x0r - x2r;
2654     y4i = x0i - x2i;
2655     y8r = x1r - x3i;
2656     y8i = x1i + x3r;
2657     y12r = x1r + x3i;
2658     y12i = x1i - x3r;
2659     x0r = a[2] + a[18];
2660     x0i = a[3] + a[19];
2661     x1r = a[2] - a[18];
2662     x1i = a[3] - a[19];
2663     x2r = a[10] + a[26];
2664     x2i = a[11] + a[27];
2665     x3r = a[10] - a[26];
2666     x3i = a[11] - a[27];
2667     y1r = x0r + x2r;
2668     y1i = x0i + x2i;
2669     y5r = x0r - x2r;
2670     y5i = x0i - x2i;
2671     x0r = x1r - x3i;
2672     x0i = x1i + x3r;
2673     y9r = wk1r * x0r - wk1i * x0i;
2674     y9i = wk1r * x0i + wk1i * x0r;
2675     x0r = x1r + x3i;
2676     x0i = x1i - x3r;
2677     y13r = wk1i * x0r - wk1r * x0i;
2678     y13i = wk1i * x0i + wk1r * x0r;
2679     x0r = a[4] + a[20];
2680     x0i = a[5] + a[21];
2681     x1r = a[4] - a[20];
2682     x1i = a[5] - a[21];
2683     x2r = a[12] + a[28];
2684     x2i = a[13] + a[29];
2685     x3r = a[12] - a[28];
2686     x3i = a[13] - a[29];
2687     y2r = x0r + x2r;
2688     y2i = x0i + x2i;
2689     y6r = x0r - x2r;
2690     y6i = x0i - x2i;
2691     x0r = x1r - x3i;
2692     x0i = x1i + x3r;
2693     y10r = wn4r * (x0r - x0i);
2694     y10i = wn4r * (x0i + x0r);
2695     x0r = x1r + x3i;
2696     x0i = x1i - x3r;
2697     y14r = wn4r * (x0r + x0i);
2698     y14i = wn4r * (x0i - x0r);
2699     x0r = a[6] + a[22];
2700     x0i = a[7] + a[23];
2701     x1r = a[6] - a[22];
2702     x1i = a[7] - a[23];
2703     x2r = a[14] + a[30];
2704     x2i = a[15] + a[31];
2705     x3r = a[14] - a[30];
2706     x3i = a[15] - a[31];
2707     y3r = x0r + x2r;
2708     y3i = x0i + x2i;
2709     y7r = x0r - x2r;
2710     y7i = x0i - x2i;
2711     x0r = x1r - x3i;
2712     x0i = x1i + x3r;
2713     y11r = wk1i * x0r - wk1r * x0i;
2714     y11i = wk1i * x0i + wk1r * x0r;
2715     x0r = x1r + x3i;
2716     x0i = x1i - x3r;
2717     y15r = wk1r * x0r - wk1i * x0i;
2718     y15i = wk1r * x0i + wk1i * x0r;
2719     x0r = y12r - y14r;
2720     x0i = y12i - y14i;
2721     x1r = y12r + y14r;
2722     x1i = y12i + y14i;
2723     x2r = y13r - y15r;
2724     x2i = y13i - y15i;
2725     x3r = y13r + y15r;
2726     x3i = y13i + y15i;
2727     a[24] = x0r + x2r;
2728     a[25] = x0i + x2i;
2729     a[26] = x0r - x2r;
2730     a[27] = x0i - x2i;
2731     a[28] = x1r - x3i;
2732     a[29] = x1i + x3r;
2733     a[30] = x1r + x3i;
2734     a[31] = x1i - x3r;
2735     x0r = y8r + y10r;
2736     x0i = y8i + y10i;
2737     x1r = y8r - y10r;
2738     x1i = y8i - y10i;
2739     x2r = y9r + y11r;
2740     x2i = y9i + y11i;
2741     x3r = y9r - y11r;
2742     x3i = y9i - y11i;
2743     a[16] = x0r + x2r;
2744     a[17] = x0i + x2i;
2745     a[18] = x0r - x2r;
2746     a[19] = x0i - x2i;
2747     a[20] = x1r - x3i;
2748     a[21] = x1i + x3r;
2749     a[22] = x1r + x3i;
2750     a[23] = x1i - x3r;
2751     x0r = y5r - y7i;
2752     x0i = y5i + y7r;
2753     x2r = wn4r * (x0r - x0i);
2754     x2i = wn4r * (x0i + x0r);
2755     x0r = y5r + y7i;
2756     x0i = y5i - y7r;
2757     x3r = wn4r * (x0r - x0i);
2758     x3i = wn4r * (x0i + x0r);
2759     x0r = y4r - y6i;
2760     x0i = y4i + y6r;
2761     x1r = y4r + y6i;
2762     x1i = y4i - y6r;
2763     a[8] = x0r + x2r;
2764     a[9] = x0i + x2i;
2765     a[10] = x0r - x2r;
2766     a[11] = x0i - x2i;
2767     a[12] = x1r - x3i;
2768     a[13] = x1i + x3r;
2769     a[14] = x1r + x3i;
2770     a[15] = x1i - x3r;
2771     x0r = y0r + y2r;
2772     x0i = y0i + y2i;
2773     x1r = y0r - y2r;
2774     x1i = y0i - y2i;
2775     x2r = y1r + y3r;
2776     x2i = y1i + y3i;
2777     x3r = y1r - y3r;
2778     x3i = y1i - y3i;
2779     a[0] = x0r + x2r;
2780     a[1] = x0i + x2i;
2781     a[2] = x0r - x2r;
2782     a[3] = x0i - x2i;
2783     a[4] = x1r - x3i;
2784     a[5] = x1i + x3r;
2785     a[6] = x1r + x3i;
2786     a[7] = x1i - x3r;
2787 }
2788 
2789 
2790 void cftf162(double *a)
2791 {
2792     double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
2793         x0r, x0i, x1r, x1i, x2r, x2i, 
2794         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
2795         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
2796         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
2797         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2798     
2799     wn4r = WR5000;
2800     wk1r = WR1250;
2801     wk1i = WI1250;
2802     wk2r = WR2500;
2803     wk2i = WI2500;
2804     wk3r = WR3750;
2805     wk3i = WI3750;
2806     x1r = a[0] - a[17];
2807     x1i = a[1] + a[16];
2808     x0r = a[8] - a[25];
2809     x0i = a[9] + a[24];
2810     x2r = wn4r * (x0r - x0i);
2811     x2i = wn4r * (x0i + x0r);
2812     y0r = x1r + x2r;
2813     y0i = x1i + x2i;
2814     y4r = x1r - x2r;
2815     y4i = x1i - x2i;
2816     x1r = a[0] + a[17];
2817     x1i = a[1] - a[16];
2818     x0r = a[8] + a[25];
2819     x0i = a[9] - a[24];
2820     x2r = wn4r * (x0r - x0i);
2821     x2i = wn4r * (x0i + x0r);
2822     y8r = x1r - x2i;
2823     y8i = x1i + x2r;
2824     y12r = x1r + x2i;
2825     y12i = x1i - x2r;
2826     x0r = a[2] - a[19];
2827     x0i = a[3] + a[18];
2828     x1r = wk1r * x0r - wk1i * x0i;
2829     x1i = wk1r * x0i + wk1i * x0r;
2830     x0r = a[10] - a[27];
2831     x0i = a[11] + a[26];
2832     x2r = wk3i * x0r - wk3r * x0i;
2833     x2i = wk3i * x0i + wk3r * x0r;
2834     y1r = x1r + x2r;
2835     y1i = x1i + x2i;
2836     y5r = x1r - x2r;
2837     y5i = x1i - x2i;
2838     x0r = a[2] + a[19];
2839     x0i = a[3] - a[18];
2840     x1r = wk3r * x0r - wk3i * x0i;
2841     x1i = wk3r * x0i + wk3i * x0r;
2842     x0r = a[10] + a[27];
2843     x0i = a[11] - a[26];
2844     x2r = wk1r * x0r + wk1i * x0i;
2845     x2i = wk1r * x0i - wk1i * x0r;
2846     y9r = x1r - x2r;
2847     y9i = x1i - x2i;
2848     y13r = x1r + x2r;
2849     y13i = x1i + x2i;
2850     x0r = a[4] - a[21];
2851     x0i = a[5] + a[20];
2852     x1r = wk2r * x0r - wk2i * x0i;
2853     x1i = wk2r * x0i + wk2i * x0r;
2854     x0r = a[12] - a[29];
2855     x0i = a[13] + a[28];
2856     x2r = wk2i * x0r - wk2r * x0i;
2857     x2i = wk2i * x0i + wk2r * x0r;
2858     y2r = x1r + x2r;
2859     y2i = x1i + x2i;
2860     y6r = x1r - x2r;
2861     y6i = x1i - x2i;
2862     x0r = a[4] + a[21];
2863     x0i = a[5] - a[20];
2864     x1r = wk2i * x0r - wk2r * x0i;
2865     x1i = wk2i * x0i + wk2r * x0r;
2866     x0r = a[12] + a[29];
2867     x0i = a[13] - a[28];
2868     x2r = wk2r * x0r - wk2i * x0i;
2869     x2i = wk2r * x0i + wk2i * x0r;
2870     y10r = x1r - x2r;
2871     y10i = x1i - x2i;
2872     y14r = x1r + x2r;
2873     y14i = x1i + x2i;
2874     x0r = a[6] - a[23];
2875     x0i = a[7] + a[22];
2876     x1r = wk3r * x0r - wk3i * x0i;
2877     x1i = wk3r * x0i + wk3i * x0r;
2878     x0r = a[14] - a[31];
2879     x0i = a[15] + a[30];
2880     x2r = wk1i * x0r - wk1r * x0i;
2881     x2i = wk1i * x0i + wk1r * x0r;
2882     y3r = x1r + x2r;
2883     y3i = x1i + x2i;
2884     y7r = x1r - x2r;
2885     y7i = x1i - x2i;
2886     x0r = a[6] + a[23];
2887     x0i = a[7] - a[22];
2888     x1r = wk1i * x0r + wk1r * x0i;
2889     x1i = wk1i * x0i - wk1r * x0r;
2890     x0r = a[14] + a[31];
2891     x0i = a[15] - a[30];
2892     x2r = wk3i * x0r - wk3r * x0i;
2893     x2i = wk3i * x0i + wk3r * x0r;
2894     y11r = x1r + x2r;
2895     y11i = x1i + x2i;
2896     y15r = x1r - x2r;
2897     y15i = x1i - x2i;
2898     x1r = y0r + y2r;
2899     x1i = y0i + y2i;
2900     x2r = y1r + y3r;
2901     x2i = y1i + y3i;
2902     a[0] = x1r + x2r;
2903     a[1] = x1i + x2i;
2904     a[2] = x1r - x2r;
2905     a[3] = x1i - x2i;
2906     x1r = y0r - y2r;
2907     x1i = y0i - y2i;
2908     x2r = y1r - y3r;
2909     x2i = y1i - y3i;
2910     a[4] = x1r - x2i;
2911     a[5] = x1i + x2r;
2912     a[6] = x1r + x2i;
2913     a[7] = x1i - x2r;
2914     x1r = y4r - y6i;
2915     x1i = y4i + y6r;
2916     x0r = y5r - y7i;
2917     x0i = y5i + y7r;
2918     x2r = wn4r * (x0r - x0i);
2919     x2i = wn4r * (x0i + x0r);
2920     a[8] = x1r + x2r;
2921     a[9] = x1i + x2i;
2922     a[10] = x1r - x2r;
2923     a[11] = x1i - x2i;
2924     x1r = y4r + y6i;
2925     x1i = y4i - y6r;
2926     x0r = y5r + y7i;
2927     x0i = y5i - y7r;
2928     x2r = wn4r * (x0r - x0i);
2929     x2i = wn4r * (x0i + x0r);
2930     a[12] = x1r - x2i;
2931     a[13] = x1i + x2r;
2932     a[14] = x1r + x2i;
2933     a[15] = x1i - x2r;
2934     x1r = y8r + y10r;
2935     x1i = y8i + y10i;
2936     x2r = y9r - y11r;
2937     x2i = y9i - y11i;
2938     a[16] = x1r + x2r;
2939     a[17] = x1i + x2i;
2940     a[18] = x1r - x2r;
2941     a[19] = x1i - x2i;
2942     x1r = y8r - y10r;
2943     x1i = y8i - y10i;
2944     x2r = y9r + y11r;
2945     x2i = y9i + y11i;
2946     a[20] = x1r - x2i;
2947     a[21] = x1i + x2r;
2948     a[22] = x1r + x2i;
2949     a[23] = x1i - x2r;
2950     x1r = y12r - y14i;
2951     x1i = y12i + y14r;
2952     x0r = y13r + y15i;
2953     x0i = y13i - y15r;
2954     x2r = wn4r * (x0r - x0i);
2955     x2i = wn4r * (x0i + x0r);
2956     a[24] = x1r + x2r;
2957     a[25] = x1i + x2i;
2958     a[26] = x1r - x2r;
2959     a[27] = x1i - x2i;
2960     x1r = y12r + y14i;
2961     x1i = y12i - y14r;
2962     x0r = y13r - y15i;
2963     x0i = y13i + y15r;
2964     x2r = wn4r * (x0r - x0i);
2965     x2i = wn4r * (x0i + x0r);
2966     a[28] = x1r - x2i;
2967     a[29] = x1i + x2r;
2968     a[30] = x1r + x2i;
2969     a[31] = x1i - x2r;
2970 }
2971 
2972 
2973 void cftf081(double *a)
2974 {
2975     double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
2976         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
2977         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
2978     
2979     wn4r = WR5000;
2980     x0r = a[0] + a[8];
2981     x0i = a[1] + a[9];
2982     x1r = a[0] - a[8];
2983     x1i = a[1] - a[9];
2984     x2r = a[4] + a[12];
2985     x2i = a[5] + a[13];
2986     x3r = a[4] - a[12];
2987     x3i = a[5] - a[13];
2988     y0r = x0r + x2r;
2989     y0i = x0i + x2i;
2990     y2r = x0r - x2r;
2991     y2i = x0i - x2i;
2992     y1r = x1r - x3i;
2993     y1i = x1i + x3r;
2994     y3r = x1r + x3i;
2995     y3i = x1i - x3r;
2996     x0r = a[2] + a[10];
2997     x0i = a[3] + a[11];
2998     x1r = a[2] - a[10];
2999     x1i = a[3] - a[11];
3000     x2r = a[6] + a[14];
3001     x2i = a[7] + a[15];
3002     x3r = a[6] - a[14];
3003     x3i = a[7] - a[15];
3004     y4r = x0r + x2r;
3005     y4i = x0i + x2i;
3006     y6r = x0r - x2r;
3007     y6i = x0i - x2i;
3008     x0r = x1r - x3i;
3009     x0i = x1i + x3r;
3010     x2r = x1r + x3i;
3011     x2i = x1i - x3r;
3012     y5r = wn4r * (x0r - x0i);
3013     y5i = wn4r * (x0r + x0i);
3014     y7r = wn4r * (x2r - x2i);
3015     y7i = wn4r * (x2r + x2i);
3016     a[8] = y1r + y5r;
3017     a[9] = y1i + y5i;
3018     a[10] = y1r - y5r;
3019     a[11] = y1i - y5i;
3020     a[12] = y3r - y7i;
3021     a[13] = y3i + y7r;
3022     a[14] = y3r + y7i;
3023     a[15] = y3i - y7r;
3024     a[0] = y0r + y4r;
3025     a[1] = y0i + y4i;
3026     a[2] = y0r - y4r;
3027     a[3] = y0i - y4i;
3028     a[4] = y2r - y6i;
3029     a[5] = y2i + y6r;
3030     a[6] = y2r + y6i;
3031     a[7] = y2i - y6r;
3032 }
3033 
3034 
3035 void cftf082(double *a)
3036 {
3037     double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
3038         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
3039         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3040     
3041     wn4r = WR5000;
3042     wk1r = WR2500;
3043     wk1i = WI2500;
3044     y0r = a[0] - a[9];
3045     y0i = a[1] + a[8];
3046     y1r = a[0] + a[9];
3047     y1i = a[1] - a[8];
3048     x0r = a[4] - a[13];
3049     x0i = a[5] + a[12];
3050     y2r = wn4r * (x0r - x0i);
3051     y2i = wn4r * (x0i + x0r);
3052     x0r = a[4] + a[13];
3053     x0i = a[5] - a[12];
3054     y3r = wn4r * (x0r - x0i);
3055     y3i = wn4r * (x0i + x0r);
3056     x0r = a[2] - a[11];
3057     x0i = a[3] + a[10];
3058     y4r = wk1r * x0r - wk1i * x0i;
3059     y4i = wk1r * x0i + wk1i * x0r;
3060     x0r = a[2] + a[11];
3061     x0i = a[3] - a[10];
3062     y5r = wk1i * x0r - wk1r * x0i;
3063     y5i = wk1i * x0i + wk1r * x0r;
3064     x0r = a[6] - a[15];
3065     x0i = a[7] + a[14];
3066     y6r = wk1i * x0r - wk1r * x0i;
3067     y6i = wk1i * x0i + wk1r * x0r;
3068     x0r = a[6] + a[15];
3069     x0i = a[7] - a[14];
3070     y7r = wk1r * x0r - wk1i * x0i;
3071     y7i = wk1r * x0i + wk1i * x0r;
3072     x0r = y0r + y2r;
3073     x0i = y0i + y2i;
3074     x1r = y4r + y6r;
3075     x1i = y4i + y6i;
3076     a[0] = x0r + x1r;
3077     a[1] = x0i + x1i;
3078     a[2] = x0r - x1r;
3079     a[3] = x0i - x1i;
3080     x0r = y0r - y2r;
3081     x0i = y0i - y2i;
3082     x1r = y4r - y6r;
3083     x1i = y4i - y6i;
3084     a[4] = x0r - x1i;
3085     a[5] = x0i + x1r;
3086     a[6] = x0r + x1i;
3087     a[7] = x0i - x1r;
3088     x0r = y1r - y3i;
3089     x0i = y1i + y3r;
3090     x1r = y5r - y7r;
3091     x1i = y5i - y7i;
3092     a[8] = x0r + x1r;
3093     a[9] = x0i + x1i;
3094     a[10] = x0r - x1r;
3095     a[11] = x0i - x1i;
3096     x0r = y1r + y3i;
3097     x0i = y1i - y3r;
3098     x1r = y5r + y7r;
3099     x1i = y5i + y7i;
3100     a[12] = x0r - x1i;
3101     a[13] = x0i + x1r;
3102     a[14] = x0r + x1i;
3103     a[15] = x0i - x1r;
3104 }
3105 
3106 
3107 void cftf040(double *a)
3108 {
3109     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3110     
3111     x0r = a[0] + a[4];
3112     x0i = a[1] + a[5];
3113     x1r = a[0] - a[4];
3114     x1i = a[1] - a[5];
3115     x2r = a[2] + a[6];
3116     x2i = a[3] + a[7];
3117     x3r = a[2] - a[6];
3118     x3i = a[3] - a[7];
3119     a[0] = x0r + x2r;
3120     a[1] = x0i + x2i;
3121     a[2] = x1r - x3i;
3122     a[3] = x1i + x3r;
3123     a[4] = x0r - x2r;
3124     a[5] = x0i - x2i;
3125     a[6] = x1r + x3i;
3126     a[7] = x1i - x3r;
3127 }
3128 
3129 
3130 void cftb040(double *a)
3131 {
3132     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3133     
3134     x0r = a[0] + a[4];
3135     x0i = a[1] + a[5];
3136     x1r = a[0] - a[4];
3137     x1i = a[1] - a[5];
3138     x2r = a[2] + a[6];
3139     x2i = a[3] + a[7];
3140     x3r = a[2] - a[6];
3141     x3i = a[3] - a[7];
3142     a[0] = x0r + x2r;
3143     a[1] = x0i + x2i;
3144     a[2] = x1r + x3i;
3145     a[3] = x1i - x3r;
3146     a[4] = x0r - x2r;
3147     a[5] = x0i - x2i;
3148     a[6] = x1r - x3i;
3149     a[7] = x1i + x3r;
3150 }
3151 
3152 
3153 void cftx020(double *a)
3154 {
3155     double x0r, x0i;
3156     
3157     x0r = a[0] - a[2];
3158     x0i = a[1] - a[3];
3159     a[0] += a[2];
3160     a[1] += a[3];
3161     a[2] = x0r;
3162     a[3] = x0i;
3163 }
3164 
3165 
3166 void rftfsub(int n, double *a)
3167 {
3168     int i, i0, j, k;
3169     double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3170     
3171     ec = 2 * M_PI_2 / n;
3172     wkr = 0;
3173     wki = 0;
3174     wdi = cos(ec);
3175     wdr = sin(ec);
3176     wdi *= wdr;
3177     wdr *= wdr;
3178     w1r = 1 - 2 * wdr;
3179     w1i = 2 * wdi;
3180     ss = 2 * w1i;
3181     i = n >> 1;
3182     for (;;) {
3183         i0 = i - 4 * RDFT_LOOP_DIV;
3184         if (i0 < 4) {
3185             i0 = 4;
3186         }
3187         for (j = i - 4; j >= i0; j -= 4) {
3188             k = n - j;
3189             xr = a[j + 2] - a[k - 2];
3190             xi = a[j + 3] + a[k - 1];
3191             yr = wdr * xr - wdi * xi;
3192             yi = wdr * xi + wdi * xr;
3193             a[j + 2] -= yr;
3194             a[j + 3] -= yi;
3195             a[k - 2] += yr;
3196             a[k - 1] -= yi;
3197             wkr += ss * wdi;
3198             wki += ss * (0.5 - wdr);
3199             xr = a[j] - a[k];
3200             xi = a[j + 1] + a[k + 1];
3201             yr = wkr * xr - wki * xi;
3202             yi = wkr * xi + wki * xr;
3203             a[j] -= yr;
3204             a[j + 1] -= yi;
3205             a[k] += yr;
3206             a[k + 1] -= yi;
3207             wdr += ss * wki;
3208             wdi += ss * (0.5 - wkr);
3209         }
3210         if (i0 == 4) {
3211             break;
3212         }
3213         wkr = 0.5 * sin(ec * i0);
3214         wki = 0.5 * cos(ec * i0);
3215         wdr = 0.5 - (wkr * w1r - wki * w1i);
3216         wdi = wkr * w1i + wki * w1r;
3217         wkr = 0.5 - wkr;
3218         i = i0;
3219     }
3220     xr = a[2] - a[n - 2];
3221     xi = a[3] + a[n - 1];
3222     yr = wdr * xr - wdi * xi;
3223     yi = wdr * xi + wdi * xr;
3224     a[2] -= yr;
3225     a[3] -= yi;
3226     a[n - 2] += yr;
3227     a[n - 1] -= yi;
3228 }
3229 
3230 
3231 void rftbsub(int n, double *a)
3232 {
3233     int i, i0, j, k;
3234     double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3235     
3236     ec = 2 * M_PI_2 / n;
3237     wkr = 0;
3238     wki = 0;
3239     wdi = cos(ec);
3240     wdr = sin(ec);
3241     wdi *= wdr;
3242     wdr *= wdr;
3243     w1r = 1 - 2 * wdr;
3244     w1i = 2 * wdi;
3245     ss = 2 * w1i;
3246     i = n >> 1;
3247     for (;;) {
3248         i0 = i - 4 * RDFT_LOOP_DIV;
3249         if (i0 < 4) {
3250             i0 = 4;
3251         }
3252         for (j = i - 4; j >= i0; j -= 4) {
3253             k = n - j;
3254             xr = a[j + 2] - a[k - 2];
3255             xi = a[j + 3] + a[k - 1];
3256             yr = wdr * xr + wdi * xi;
3257             yi = wdr * xi - wdi * xr;
3258             a[j + 2] -= yr;
3259             a[j + 3] -= yi;
3260             a[k - 2] += yr;
3261             a[k - 1] -= yi;
3262             wkr += ss * wdi;
3263             wki += ss * (0.5 - wdr);
3264             xr = a[j] - a[k];
3265             xi = a[j + 1] + a[k + 1];
3266             yr = wkr * xr + wki * xi;
3267             yi = wkr * xi - wki * xr;
3268             a[j] -= yr;
3269             a[j + 1] -= yi;
3270             a[k] += yr;
3271             a[k + 1] -= yi;
3272             wdr += ss * wki;
3273             wdi += ss * (0.5 - wkr);
3274         }
3275         if (i0 == 4) {
3276             break;
3277         }
3278         wkr = 0.5 * sin(ec * i0);
3279         wki = 0.5 * cos(ec * i0);
3280         wdr = 0.5 - (wkr * w1r - wki * w1i);
3281         wdi = wkr * w1i + wki * w1r;
3282         wkr = 0.5 - wkr;
3283         i = i0;
3284     }
3285     xr = a[2] - a[n - 2];
3286     xi = a[3] + a[n - 1];
3287     yr = wdr * xr + wdi * xi;
3288     yi = wdr * xi - wdi * xr;
3289     a[2] -= yr;
3290     a[3] -= yi;
3291     a[n - 2] += yr;
3292     a[n - 1] -= yi;
3293 }
3294 
3295 
3296 void dctsub(int n, double *a)
3297 {
3298     int i, i0, j, k, m;
3299     double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3300     
3301     ec = M_PI_2 / n;
3302     wkr = 0.5;
3303     wki = 0.5;
3304     w1r = cos(ec);
3305     w1i = sin(ec);
3306     wdr = 0.5 * (w1r - w1i);
3307     wdi = 0.5 * (w1r + w1i);
3308     ss = 2 * w1i;
3309     m = n >> 1;
3310     i = 0;
3311     for (;;) {
3312         i0 = i + 2 * DCST_LOOP_DIV;
3313         if (i0 > m - 2) {
3314             i0 = m - 2;
3315         }
3316         for (j = i + 2; j <= i0; j += 2) {
3317             k = n - j;
3318             xr = wdi * a[j - 1] - wdr * a[k + 1];
3319             xi = wdr * a[j - 1] + wdi * a[k + 1];
3320             wkr -= ss * wdi;
3321             wki += ss * wdr;
3322             yr = wki * a[j] - wkr * a[k];
3323             yi = wkr * a[j] + wki * a[k];
3324             wdr -= ss * wki;
3325             wdi += ss * wkr;
3326             a[k + 1] = xr;
3327             a[k] = yr;
3328             a[j - 1] = xi;
3329             a[j] = yi;
3330         }
3331         if (i0 == m - 2) {
3332             break;
3333         }
3334         wdr = cos(ec * i0);
3335         wdi = sin(ec * i0);
3336         wkr = 0.5 * (wdr - wdi);
3337         wki = 0.5 * (wdr + wdi);
3338         wdr = wkr * w1r - wki * w1i;
3339         wdi = wkr * w1i + wki * w1r;
3340         i = i0;
3341     }
3342     xr = wdi * a[m - 1] - wdr * a[m + 1];
3343     a[m - 1] = wdr * a[m - 1] + wdi * a[m + 1];
3344     a[m + 1] = xr;
3345     a[m] *= WR5000;
3346 }
3347 
3348 
3349 void dstsub(int n, double *a)
3350 {
3351     int i, i0, j, k, m;
3352     double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3353     
3354     ec = M_PI_2 / n;
3355     wkr = 0.5;
3356     wki = 0.5;
3357     w1r = cos(ec);
3358     w1i = sin(ec);
3359     wdr = 0.5 * (w1r - w1i);
3360     wdi = 0.5 * (w1r + w1i);
3361     ss = 2 * w1i;
3362     m = n >> 1;
3363     i = 0;
3364     for (;;) {
3365         i0 = i + 2 * DCST_LOOP_DIV;
3366         if (i0 > m - 2) {
3367             i0 = m - 2;
3368         }
3369         for (j = i + 2; j <= i0; j += 2) {
3370             k = n - j;
3371             xr = wdi * a[k + 1] - wdr * a[j - 1];
3372             xi = wdr * a[k + 1] + wdi * a[j - 1];
3373             wkr -= ss * wdi;
3374             wki += ss * wdr;
3375             yr = wki * a[k] - wkr * a[j];
3376             yi = wkr * a[k] + wki * a[j];
3377             wdr -= ss * wki;
3378             wdi += ss * wkr;
3379             a[j - 1] = xr;
3380             a[j] = yr;
3381             a[k + 1] = xi;
3382             a[k] = yi;
3383         }
3384         if (i0 == m - 2) {
3385             break;
3386         }
3387         wdr = cos(ec * i0);
3388         wdi = sin(ec * i0);
3389         wkr = 0.5 * (wdr - wdi);
3390         wki = 0.5 * (wdr + wdi);
3391         wdr = wkr * w1r - wki * w1i;
3392         wdi = wkr * w1i + wki * w1r;
3393         i = i0;
3394     }
3395     xr = wdi * a[m + 1] - wdr * a[m - 1];
3396     a[m + 1] = wdr * a[m + 1] + wdi * a[m - 1];
3397     a[m - 1] = xr;
3398     a[m] *= WR5000;
3399 }
3400 
3401 
3402 void dctsub4(int n, double *a)
3403 {
3404     int m;
3405     double wki, wdr, wdi, xr;
3406     
3407     wki = WR5000;
3408     m = n >> 1;
3409     if (m == 2) {
3410         wdr = wki * WI2500;
3411         wdi = wki * WR2500;
3412         xr = wdi * a[1] - wdr * a[3];
3413         a[1] = wdr * a[1] + wdi * a[3];
3414         a[3] = xr;
3415     }
3416     a[m] *= wki;
3417 }
3418 
3419 
3420 void dstsub4(int n, double *a)
3421 {
3422     int m;
3423     double wki, wdr, wdi, xr;
3424     
3425     wki = WR5000;
3426     m = n >> 1;
3427     if (m == 2) {
3428         wdr = wki * WI2500;
3429         wdi = wki * WR2500;
3430         xr = wdi * a[3] - wdr * a[1];
3431         a[3] = wdr * a[3] + wdi * a[1];
3432         a[1] = xr;
3433     }
3434     a[m] *= wki;
3435 }
3436