home *** CD-ROM | disk | FTP | other *** search
/ Cubase Magazine 30 / Issue #30.iso / pc / 2-SOFTWARE / RTEQ 3.0 / SRC / FFTSG.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1999-06-13  |  68.9 KB  |  2,643 lines
  1. #include "precomp.h"
  2. /*
  3. Fast Fourier/Cosine/Sine Transform
  4.     dimension   :one
  5.     data length :power of 2
  6.     decimation  :frequency
  7.     radix       :split-radix
  8.     data        :inplace
  9.     table       :use
  10. functions
  11.     cdft: Complex Discrete Fourier Transform
  12.     rdft: Real Discrete Fourier Transform
  13.     ddct: Discrete Cosine Transform
  14.     ddst: Discrete Sine Transform
  15.     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
  16.     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
  17. function prototypes
  18.     void cdft(int, int, double *, int *, double *);
  19.     void rdft(int, int, double *, int *, double *);
  20.     void ddct(int, int, double *, int *, double *);
  21.     void ddst(int, int, double *, int *, double *);
  22.     void dfct(int, double *, double *, int *, double *);
  23.     void dfst(int, double *, double *, int *, double *);
  24.  
  25.  
  26. -------- Complex DFT (Discrete Fourier Transform) --------
  27.     [definition]
  28.         <case1>
  29.             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
  30.         <case2>
  31.             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
  32.         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
  33.     [usage]
  34.         <case1>
  35.             ip[0] = 0; // first time only
  36.             cdft(2*n, 1, a, ip, w);
  37.         <case2>
  38.             ip[0] = 0; // first time only
  39.             cdft(2*n, -1, a, ip, w);
  40.     [parameters]
  41.         2*n            :data length (int)
  42.                         n >= 1, n = power of 2
  43.         a[0...2*n-1]   :input/output data (double *)
  44.                         input data
  45.                             a[2*j] = Re(x[j]), 
  46.                             a[2*j+1] = Im(x[j]), 0<=j<n
  47.                         output data
  48.                             a[2*k] = Re(X[k]), 
  49.                             a[2*k+1] = Im(X[k]), 0<=k<n
  50.         ip[0...*]      :work area for bit reversal (int *)
  51.                         length of ip >= 2+sqrt(n)
  52.                         strictly, 
  53.                         length of ip >= 
  54.                             2+(1<<(int)(log(n+0.5)/log(2))/2).
  55.                         ip[0],ip[1] are pointers of the cos/sin table.
  56.         w[0...n/2-1]   :cos/sin table (double *)
  57.                         w[],ip[] are initialized if ip[0] == 0.
  58.     [remark]
  59.         Inverse of 
  60.             cdft(2*n, -1, a, ip, w);
  61.         is 
  62.             cdft(2*n, 1, a, ip, w);
  63.             for (j = 0; j <= 2 * n - 1; j++) {
  64.                 a[j] *= 1.0 / n;
  65.             }
  66.         .
  67.  
  68.  
  69. -------- Real DFT / Inverse of Real DFT --------
  70.     [definition]
  71.         <case1> RDFT
  72.             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
  73.             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
  74.         <case2> IRDFT (excluding scale)
  75.             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
  76.                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
  77.                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
  78.     [usage]
  79.         <case1>
  80.             ip[0] = 0; // first time only
  81.             rdft(n, 1, a, ip, w);
  82.         <case2>
  83.             ip[0] = 0; // first time only
  84.             rdft(n, -1, a, ip, w);
  85.     [parameters]
  86.         n              :data length (int)
  87.                         n >= 2, n = power of 2
  88.         a[0...n-1]     :input/output data (double *)
  89.                         <case1>
  90.                             output data
  91.                                 a[2*k] = R[k], 0<=k<n/2
  92.                                 a[2*k+1] = I[k], 0<k<n/2
  93.                                 a[1] = R[n/2]
  94.                         <case2>
  95.                             input data
  96.                                 a[2*j] = R[j], 0<=j<n/2
  97.                                 a[2*j+1] = I[j], 0<j<n/2
  98.                                 a[1] = R[n/2]
  99.         ip[0...*]      :work area for bit reversal (int *)
  100.                         length of ip >= 2+sqrt(n/2)
  101.                         strictly, 
  102.                         length of ip >= 
  103.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  104.                         ip[0],ip[1] are pointers of the cos/sin table.
  105.         w[0...n/2-1]   :cos/sin table (double *)
  106.                         w[],ip[] are initialized if ip[0] == 0.
  107.     [remark]
  108.         Inverse of 
  109.             rdft(n, 1, a, ip, w);
  110.         is 
  111.             rdft(n, -1, a, ip, w);
  112.             for (j = 0; j <= n - 1; j++) {
  113.                 a[j] *= 2.0 / n;
  114.             }
  115.         .
  116.  
  117.  
  118. -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
  119.     [definition]
  120.         <case1> IDCT (excluding scale)
  121.             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
  122.         <case2> DCT
  123.             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
  124.     [usage]
  125.         <case1>
  126.             ip[0] = 0; // first time only
  127.             ddct(n, 1, a, ip, w);
  128.         <case2>
  129.             ip[0] = 0; // first time only
  130.             ddct(n, -1, a, ip, w);
  131.     [parameters]
  132.         n              :data length (int)
  133.                         n >= 2, n = power of 2
  134.         a[0...n-1]     :input/output data (double *)
  135.                         output data
  136.                             a[k] = C[k], 0<=k<n
  137.         ip[0...*]      :work area for bit reversal (int *)
  138.                         length of ip >= 2+sqrt(n/2)
  139.                         strictly, 
  140.                         length of ip >= 
  141.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  142.                         ip[0],ip[1] are pointers of the cos/sin table.
  143.         w[0...n*5/4-1] :cos/sin table (double *)
  144.                         w[],ip[] are initialized if ip[0] == 0.
  145.     [remark]
  146.         Inverse of 
  147.             ddct(n, -1, a, ip, w);
  148.         is 
  149.             a[0] *= 0.5;
  150.             ddct(n, 1, a, ip, w);
  151.             for (j = 0; j <= n - 1; j++) {
  152.                 a[j] *= 2.0 / n;
  153.             }
  154.         .
  155.  
  156.  
  157. -------- DST (Discrete Sine Transform) / Inverse of DST --------
  158.     [definition]
  159.         <case1> IDST (excluding scale)
  160.             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
  161.         <case2> DST
  162.             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
  163.     [usage]
  164.         <case1>
  165.             ip[0] = 0; // first time only
  166.             ddst(n, 1, a, ip, w);
  167.         <case2>
  168.             ip[0] = 0; // first time only
  169.             ddst(n, -1, a, ip, w);
  170.     [parameters]
  171.         n              :data length (int)
  172.                         n >= 2, n = power of 2
  173.         a[0...n-1]     :input/output data (double *)
  174.                         <case1>
  175.                             input data
  176.                                 a[j] = A[j], 0<j<n
  177.                                 a[0] = A[n]
  178.                             output data
  179.                                 a[k] = S[k], 0<=k<n
  180.                         <case2>
  181.                             output data
  182.                                 a[k] = S[k], 0<k<n
  183.                                 a[0] = S[n]
  184.         ip[0...*]      :work area for bit reversal (int *)
  185.                         length of ip >= 2+sqrt(n/2)
  186.                         strictly, 
  187.                         length of ip >= 
  188.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  189.                         ip[0],ip[1] are pointers of the cos/sin table.
  190.         w[0...n*5/4-1] :cos/sin table (double *)
  191.                         w[],ip[] are initialized if ip[0] == 0.
  192.     [remark]
  193.         Inverse of 
  194.             ddst(n, -1, a, ip, w);
  195.         is 
  196.             a[0] *= 0.5;
  197.             ddst(n, 1, a, ip, w);
  198.             for (j = 0; j <= n - 1; j++) {
  199.                 a[j] *= 2.0 / n;
  200.             }
  201.         .
  202.  
  203.  
  204. -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
  205.     [definition]
  206.         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
  207.     [usage]
  208.         ip[0] = 0; // first time only
  209.         dfct(n, a, t, ip, w);
  210.     [parameters]
  211.         n              :data length - 1 (int)
  212.                         n >= 2, n = power of 2
  213.         a[0...n]       :input/output data (double *)
  214.                         output data
  215.                             a[k] = C[k], 0<=k<=n
  216.         t[0...n/2]     :work area (double *)
  217.         ip[0...*]      :work area for bit reversal (int *)
  218.                         length of ip >= 2+sqrt(n/4)
  219.                         strictly, 
  220.                         length of ip >= 
  221.                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
  222.                         ip[0],ip[1] are pointers of the cos/sin table.
  223.         w[0...n*5/8-1] :cos/sin table (double *)
  224.                         w[],ip[] are initialized if ip[0] == 0.
  225.     [remark]
  226.         Inverse of 
  227.             a[0] *= 0.5;
  228.             a[n] *= 0.5;
  229.             dfct(n, a, t, ip, w);
  230.         is 
  231.             a[0] *= 0.5;
  232.             a[n] *= 0.5;
  233.             dfct(n, a, t, ip, w);
  234.             for (j = 0; j <= n; j++) {
  235.                 a[j] *= 2.0 / n;
  236.             }
  237.         .
  238.  
  239.  
  240. -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
  241.     [definition]
  242.         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
  243.     [usage]
  244.         ip[0] = 0; // first time only
  245.         dfst(n, a, t, ip, w);
  246.     [parameters]
  247.         n              :data length + 1 (int)
  248.                         n >= 2, n = power of 2
  249.         a[0...n-1]     :input/output data (double *)
  250.                         output data
  251.                             a[k] = S[k], 0<k<n
  252.                         (a[0] is used for work area)
  253.         t[0...n/2-1]   :work area (double *)
  254.         ip[0...*]      :work area for bit reversal (int *)
  255.                         length of ip >= 2+sqrt(n/4)
  256.                         strictly, 
  257.                         length of ip >= 
  258.                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
  259.                         ip[0],ip[1] are pointers of the cos/sin table.
  260.         w[0...n*5/8-1] :cos/sin table (double *)
  261.                         w[],ip[] are initialized if ip[0] == 0.
  262.     [remark]
  263.         Inverse of 
  264.             dfst(n, a, t, ip, w);
  265.         is 
  266.             dfst(n, a, t, ip, w);
  267.             for (j = 1; j <= n - 1; j++) {
  268.                 a[j] *= 2.0 / n;
  269.             }
  270.         .
  271.  
  272.  
  273. Appendix :
  274.     The cos/sin table is recalculated when the larger table required.
  275.     w[] and ip[] are compatible with all routines.
  276. */
  277.  
  278.  
  279. void cdft(int n, int isgn, double *a, int *ip, double *w)
  280. {
  281.     void makewt(int nw, int *ip, double *w);
  282.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  283.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  284.     int nw;
  285.     
  286.     nw = ip[0];
  287.     if (n > (nw << 2)) {
  288.         nw = n >> 2;
  289.         makewt(nw, ip, w);
  290.     }
  291.     if (isgn >= 0) {
  292.         cftfsub(n, a, ip + 2, nw, w);
  293.     } else {
  294.         cftbsub(n, a, ip + 2, nw, w);
  295.     }
  296. }
  297.  
  298.  
  299. void rdft(int n, int isgn, double *a, int *ip, double *w)
  300. {
  301.     void makewt(int nw, int *ip, double *w);
  302.     void makect(int nc, int *ip, double *c);
  303.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  304.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  305.     void rftfsub(int n, double *a, int nc, double *c);
  306.     void rftbsub(int n, double *a, int nc, double *c);
  307.     int nw, nc;
  308.     double xi;
  309.     
  310.     nw = ip[0];
  311.     if (n > (nw << 2)) {
  312.         nw = n >> 2;
  313.         makewt(nw, ip, w);
  314.     }
  315.     nc = ip[1];
  316.     if (n > (nc << 2)) {
  317.         nc = n >> 2;
  318.         makect(nc, ip, w + nw);
  319.     }
  320.     if (isgn >= 0) {
  321.         if (n > 4) {
  322.             cftfsub(n, a, ip + 2, nw, w);
  323.             rftfsub(n, a, nc, w + nw);
  324.         } else if (n == 4) {
  325.             cftfsub(n, a, ip + 2, nw, w);
  326.         }
  327.         xi = a[0] - a[1];
  328.         a[0] += a[1];
  329.         a[1] = xi;
  330.     } else {
  331.         a[1] = 0.5 * (a[0] - a[1]);
  332.         a[0] -= a[1];
  333.         if (n > 4) {
  334.             rftbsub(n, a, nc, w + nw);
  335.             cftbsub(n, a, ip + 2, nw, w);
  336.         } else if (n == 4) {
  337.             cftbsub(n, a, ip + 2, nw, w);
  338.         }
  339.     }
  340. }
  341.  
  342.  
  343. void ddct(int n, int isgn, double *a, int *ip, double *w)
  344. {
  345.     void makewt(int nw, int *ip, double *w);
  346.     void makect(int nc, int *ip, double *c);
  347.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  348.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  349.     void rftfsub(int n, double *a, int nc, double *c);
  350.     void rftbsub(int n, double *a, int nc, double *c);
  351.     void dctsub(int n, double *a, int nc, double *c);
  352.     int j, nw, nc;
  353.     double xr;
  354.     
  355.     nw = ip[0];
  356.     if (n > (nw << 2)) {
  357.         nw = n >> 2;
  358.         makewt(nw, ip, w);
  359.     }
  360.     nc = ip[1];
  361.     if (n > nc) {
  362.         nc = n;
  363.         makect(nc, ip, w + nw);
  364.     }
  365.     if (isgn < 0) {
  366.         xr = a[n - 1];
  367.         for (j = n - 2; j >= 2; j -= 2) {
  368.             a[j + 1] = a[j] - a[j - 1];
  369.             a[j] += a[j - 1];
  370.         }
  371.         a[1] = a[0] - xr;
  372.         a[0] += xr;
  373.         if (n > 4) {
  374.             rftbsub(n, a, nc, w + nw);
  375.             cftbsub(n, a, ip + 2, nw, w);
  376.         } else if (n == 4) {
  377.             cftbsub(n, a, ip + 2, nw, w);
  378.         }
  379.     }
  380.     dctsub(n, a, nc, w + nw);
  381.     if (isgn >= 0) {
  382.         if (n > 4) {
  383.             cftfsub(n, a, ip + 2, nw, w);
  384.             rftfsub(n, a, nc, w + nw);
  385.         } else if (n == 4) {
  386.             cftfsub(n, a, ip + 2, nw, w);
  387.         }
  388.         xr = a[0] - a[1];
  389.         a[0] += a[1];
  390.         for (j = 2; j < n; j += 2) {
  391.             a[j - 1] = a[j] - a[j + 1];
  392.             a[j] += a[j + 1];
  393.         }
  394.         a[n - 1] = xr;
  395.     }
  396. }
  397.  
  398.  
  399. void ddst(int n, int isgn, double *a, int *ip, double *w)
  400. {
  401.     void makewt(int nw, int *ip, double *w);
  402.     void makect(int nc, int *ip, double *c);
  403.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  404.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  405.     void rftfsub(int n, double *a, int nc, double *c);
  406.     void rftbsub(int n, double *a, int nc, double *c);
  407.     void dstsub(int n, double *a, int nc, double *c);
  408.     int j, nw, nc;
  409.     double xr;
  410.     
  411.     nw = ip[0];
  412.     if (n > (nw << 2)) {
  413.         nw = n >> 2;
  414.         makewt(nw, ip, w);
  415.     }
  416.     nc = ip[1];
  417.     if (n > nc) {
  418.         nc = n;
  419.         makect(nc, ip, w + nw);
  420.     }
  421.     if (isgn < 0) {
  422.         xr = a[n - 1];
  423.         for (j = n - 2; j >= 2; j -= 2) {
  424.             a[j + 1] = -a[j] - a[j - 1];
  425.             a[j] -= a[j - 1];
  426.         }
  427.         a[1] = a[0] + xr;
  428.         a[0] -= xr;
  429.         if (n > 4) {
  430.             rftbsub(n, a, nc, w + nw);
  431.             cftbsub(n, a, ip + 2, nw, w);
  432.         } else if (n == 4) {
  433.             cftbsub(n, a, ip + 2, nw, w);
  434.         }
  435.     }
  436.     dstsub(n, a, nc, w + nw);
  437.     if (isgn >= 0) {
  438.         if (n > 4) {
  439.             cftfsub(n, a, ip + 2, nw, w);
  440.             rftfsub(n, a, nc, w + nw);
  441.         } else if (n == 4) {
  442.             cftfsub(n, a, ip + 2, nw, w);
  443.         }
  444.         xr = a[0] - a[1];
  445.         a[0] += a[1];
  446.         for (j = 2; j < n; j += 2) {
  447.             a[j - 1] = -a[j] - a[j + 1];
  448.             a[j] -= a[j + 1];
  449.         }
  450.         a[n - 1] = -xr;
  451.     }
  452. }
  453.  
  454.  
  455. void dfct(int n, double *a, double *t, int *ip, double *w)
  456. {
  457.     void makewt(int nw, int *ip, double *w);
  458.     void makect(int nc, int *ip, double *c);
  459.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  460.     void rftfsub(int n, double *a, int nc, double *c);
  461.     void dctsub(int n, double *a, int nc, double *c);
  462.     int j, k, l, m, mh, nw, nc;
  463.     double xr, xi, yr, yi;
  464.     
  465.     nw = ip[0];
  466.     if (n > (nw << 3)) {
  467.         nw = n >> 3;
  468.         makewt(nw, ip, w);
  469.     }
  470.     nc = ip[1];
  471.     if (n > (nc << 1)) {
  472.         nc = n >> 1;
  473.         makect(nc, ip, w + nw);
  474.     }
  475.     m = n >> 1;
  476.     yi = a[m];
  477.     xi = a[0] + a[n];
  478.     a[0] -= a[n];
  479.     t[0] = xi - yi;
  480.     t[m] = xi + yi;
  481.     if (n > 2) {
  482.         mh = m >> 1;
  483.         for (j = 1; j < mh; j++) {
  484.             k = m - j;
  485.             xr = a[j] - a[n - j];
  486.             xi = a[j] + a[n - j];
  487.             yr = a[k] - a[n - k];
  488.             yi = a[k] + a[n - k];
  489.             a[j] = xr;
  490.             a[k] = yr;
  491.             t[j] = xi - yi;
  492.             t[k] = xi + yi;
  493.         }
  494.         t[mh] = a[mh] + a[n - mh];
  495.         a[mh] -= a[n - mh];
  496.         dctsub(m, a, nc, w + nw);
  497.         if (m > 4) {
  498.             cftfsub(m, a, ip + 2, nw, w);
  499.             rftfsub(m, a, nc, w + nw);
  500.         } else if (m == 4) {
  501.             cftfsub(m, a, ip + 2, nw, w);
  502.         }
  503.         a[n - 1] = a[0] - a[1];
  504.         a[1] = a[0] + a[1];
  505.         for (j = m - 2; j >= 2; j -= 2) {
  506.             a[2 * j + 1] = a[j] + a[j + 1];
  507.             a[2 * j - 1] = a[j] - a[j + 1];
  508.         }
  509.         l = 2;
  510.         m = mh;
  511.         while (m >= 2) {
  512.             dctsub(m, t, nc, w + nw);
  513.             if (m > 4) {
  514.                 cftfsub(m, t, ip + 2, nw, w);
  515.                 rftfsub(m, t, nc, w + nw);
  516.             } else if (m == 4) {
  517.                 cftfsub(m, t, ip + 2, nw, w);
  518.             }
  519.             a[n - l] = t[0] - t[1];
  520.             a[l] = t[0] + t[1];
  521.             k = 0;
  522.             for (j = 2; j < m; j += 2) {
  523.                 k += l << 2;
  524.                 a[k - l] = t[j] - t[j + 1];
  525.                 a[k + l] = t[j] + t[j + 1];
  526.             }
  527.             l <<= 1;
  528.             mh = m >> 1;
  529.             for (j = 0; j < mh; j++) {
  530.                 k = m - j;
  531.                 t[j] = t[m + k] - t[m + j];
  532.                 t[k] = t[m + k] + t[m + j];
  533.             }
  534.             t[mh] = t[m + mh];
  535.             m = mh;
  536.         }
  537.         a[l] = t[0];
  538.         a[n] = t[2] - t[1];
  539.         a[0] = t[2] + t[1];
  540.     } else {
  541.         a[1] = a[0];
  542.         a[2] = t[0];
  543.         a[0] = t[1];
  544.     }
  545. }
  546.  
  547.  
  548. void dfst(int n, double *a, double *t, int *ip, double *w)
  549. {
  550.     void makewt(int nw, int *ip, double *w);
  551.     void makect(int nc, int *ip, double *c);
  552.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  553.     void rftfsub(int n, double *a, int nc, double *c);
  554.     void dstsub(int n, double *a, int nc, double *c);
  555.     int j, k, l, m, mh, nw, nc;
  556.     double xr, xi, yr, yi;
  557.     
  558.     nw = ip[0];
  559.     if (n > (nw << 3)) {
  560.         nw = n >> 3;
  561.         makewt(nw, ip, w);
  562.     }
  563.     nc = ip[1];
  564.     if (n > (nc << 1)) {
  565.         nc = n >> 1;
  566.         makect(nc, ip, w + nw);
  567.     }
  568.     if (n > 2) {
  569.         m = n >> 1;
  570.         mh = m >> 1;
  571.         for (j = 1; j < mh; j++) {
  572.             k = m - j;
  573.             xr = a[j] + a[n - j];
  574.             xi = a[j] - a[n - j];
  575.             yr = a[k] + a[n - k];
  576.             yi = a[k] - a[n - k];
  577.             a[j] = xr;
  578.             a[k] = yr;
  579.             t[j] = xi + yi;
  580.             t[k] = xi - yi;
  581.         }
  582.         t[0] = a[mh] - a[n - mh];
  583.         a[mh] += a[n - mh];
  584.         a[0] = a[m];
  585.         dstsub(m, a, nc, w + nw);
  586.         if (m > 4) {
  587.             cftfsub(m, a, ip + 2, nw, w);
  588.             rftfsub(m, a, nc, w + nw);
  589.         } else if (m == 4) {
  590.             cftfsub(m, a, ip + 2, nw, w);
  591.         }
  592.         a[n - 1] = a[1] - a[0];
  593.         a[1] = a[0] + a[1];
  594.         for (j = m - 2; j >= 2; j -= 2) {
  595.             a[2 * j + 1] = a[j] - a[j + 1];
  596.             a[2 * j - 1] = -a[j] - a[j + 1];
  597.         }
  598.         l = 2;
  599.         m = mh;
  600.         while (m >= 2) {
  601.             dstsub(m, t, nc, w + nw);
  602.             if (m > 4) {
  603.                 cftfsub(m, t, ip + 2, nw, w);
  604.                 rftfsub(m, t, nc, w + nw);
  605.             } else if (m == 4) {
  606.                 cftfsub(m, t, ip + 2, nw, w);
  607.             }
  608.             a[n - l] = t[1] - t[0];
  609.             a[l] = t[0] + t[1];
  610.             k = 0;
  611.             for (j = 2; j < m; j += 2) {
  612.                 k += l << 2;
  613.                 a[k - l] = -t[j] - t[j + 1];
  614.                 a[k + l] = t[j] - t[j + 1];
  615.             }
  616.             l <<= 1;
  617.             mh = m >> 1;
  618.             for (j = 1; j < mh; j++) {
  619.                 k = m - j;
  620.                 t[j] = t[m + k] + t[m + j];
  621.                 t[k] = t[m + k] - t[m + j];
  622.             }
  623.             t[0] = t[m + mh];
  624.             m = mh;
  625.         }
  626.         a[l] = t[0];
  627.     }
  628.     a[0] = 0;
  629. }
  630.  
  631.  
  632. /* -------- initializing routines -------- */
  633.  
  634.  
  635. #include <math.h>
  636.  
  637. void makewt(int nw, int *ip, double *w)
  638. {
  639.     int j, nwh, nw0, nw1;
  640.     double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
  641.     
  642.     ip[0] = nw;
  643.     ip[1] = 1;
  644.     if (nw > 2) {
  645.         nwh = nw >> 1;
  646.         delta = atan(1.0) / nwh;
  647.         wn4r = cos(delta * nwh);
  648.         w[0] = 1;
  649.         w[1] = wn4r;
  650.         if (nwh >= 4) {
  651.             w[2] = 0.5 / cos(delta * 2);
  652.             w[3] = 0.5 / cos(delta * 6);
  653.         }
  654.         for (j = 4; j < nwh; j += 4) {
  655.             w[j] = cos(delta * j);
  656.             w[j + 1] = sin(delta * j);
  657.             w[j + 2] = cos(3 * delta * j);
  658.             w[j + 3] = sin(3 * delta * j);
  659.         }
  660.         nw0 = 0;
  661.         while (nwh > 2) {
  662.             nw1 = nw0 + nwh;
  663.             nwh >>= 1;
  664.             w[nw1] = 1;
  665.             w[nw1 + 1] = wn4r;
  666.             if (nwh >= 4) {
  667.                 wk1r = w[nw0 + 4];
  668.                 wk3r = w[nw0 + 6];
  669.                 w[nw1 + 2] = 0.5 / wk1r;
  670.                 w[nw1 + 3] = 0.5 / wk3r;
  671.             }
  672.             for (j = 4; j < nwh; j += 4) {
  673.                 wk1r = w[nw0 + 2 * j];
  674.                 wk1i = w[nw0 + 2 * j + 1];
  675.                 wk3r = w[nw0 + 2 * j + 2];
  676.                 wk3i = w[nw0 + 2 * j + 3];
  677.                 w[nw1 + j] = wk1r;
  678.                 w[nw1 + j + 1] = wk1i;
  679.                 w[nw1 + j + 2] = wk3r;
  680.                 w[nw1 + j + 3] = wk3i;
  681.             }
  682.             nw0 = nw1;
  683.         }
  684.     }
  685. }
  686.  
  687.  
  688. void makect(int nc, int *ip, double *c)
  689. {
  690.     int j, nch;
  691.     double delta;
  692.     
  693.     ip[1] = nc;
  694.     if (nc > 1) {
  695.         nch = nc >> 1;
  696.         delta = atan(1.0) / nch;
  697.         c[0] = cos(delta * nch);
  698.         c[nch] = 0.5 * c[0];
  699.         for (j = 1; j < nch; j++) {
  700.             c[j] = 0.5 * cos(delta * j);
  701.             c[nc - j] = 0.5 * sin(delta * j);
  702.         }
  703.     }
  704. }
  705.  
  706.  
  707. /* -------- child routines -------- */
  708.  
  709.  
  710. #ifndef CDFT_RECURSIVE_N  /* length of the recursive FFT mode */
  711. #define CDFT_RECURSIVE_N 512  /* <= (L1 cache size) / 16 */
  712. #endif
  713.  
  714.  
  715. void cftfsub(int n, double *a, int *ip, int nw, double *w)
  716. {
  717.     void bitrv2(int n, int *ip, double *a);
  718.     void bitrv216(double *a);
  719.     void bitrv208(double *a);
  720.     void cftf1st(int n, double *a, double *w);
  721.     void cftrec1(int n, double *a, int nw, double *w);
  722.     void cftrec2(int n, double *a, int nw, double *w);
  723.     void cftexp1(int n, double *a, int nw, double *w);
  724.     void cftfx41(int n, double *a, int nw, double *w);
  725.     void cftf161(double *a, double *w);
  726.     void cftf081(double *a, double *w);
  727.     void cftf040(double *a);
  728.     void cftx020(double *a);
  729.     int m;
  730.     
  731.     if (n > 32) {
  732.         m = n >> 2;
  733.         cftf1st(n, a, &w[nw - m]);
  734.         if (n > CDFT_RECURSIVE_N) {
  735.             cftrec1(m, a, nw, w);
  736.             cftrec2(m, &a[m], nw, w);
  737.             cftrec1(m, &a[2 * m], nw, w);
  738.             cftrec1(m, &a[3 * m], nw, w);
  739.         } else if (m > 32) {
  740.             cftexp1(n, a, nw, w);
  741.         } else {
  742.             cftfx41(n, a, nw, w);
  743.         }
  744.         bitrv2(n, ip, a);
  745.     } else if (n > 8) {
  746.         if (n == 32) {
  747.             cftf161(a, &w[nw - 8]);
  748.             bitrv216(a);
  749.         } else {
  750.             cftf081(a, w);
  751.             bitrv208(a);
  752.         }
  753.     } else if (n == 8) {
  754.         cftf040(a);
  755.     } else if (n == 4) {
  756.         cftx020(a);
  757.     }
  758. }
  759.  
  760.  
  761. void cftbsub(int n, double *a, int *ip, int nw, double *w)
  762. {
  763.     void bitrv2conj(int n, int *ip, double *a);
  764.     void bitrv216neg(double *a);
  765.     void bitrv208neg(double *a);
  766.     void cftb1st(int n, double *a, double *w);
  767.     void cftrec1(int n, double *a, int nw, double *w);
  768.     void cftrec2(int n, double *a, int nw, double *w);
  769.     void cftexp1(int n, double *a, int nw, double *w);
  770.     void cftfx41(int n, double *a, int nw, double *w);
  771.     void cftf161(double *a, double *w);
  772.     void cftf081(double *a, double *w);
  773.     void cftb040(double *a);
  774.     void cftx020(double *a);
  775.     int m;
  776.     
  777.     if (n > 32) {
  778.         m = n >> 2;
  779.         cftb1st(n, a, &w[nw - m]);
  780.         if (n > CDFT_RECURSIVE_N) {
  781.             cftrec1(m, a, nw, w);
  782.             cftrec2(m, &a[m], nw, w);
  783.             cftrec1(m, &a[2 * m], nw, w);
  784.             cftrec1(m, &a[3 * m], nw, w);
  785.         } else if (m > 32) {
  786.             cftexp1(n, a, nw, w);
  787.         } else {
  788.             cftfx41(n, a, nw, w);
  789.         }
  790.         bitrv2conj(n, ip, a);
  791.     } else if (n > 8) {
  792.         if (n == 32) {
  793.             cftf161(a, &w[nw - 8]);
  794.             bitrv216neg(a);
  795.         } else {
  796.             cftf081(a, w);
  797.             bitrv208neg(a);
  798.         }
  799.     } else if (n == 8) {
  800.         cftb040(a);
  801.     } else if (n == 4) {
  802.         cftx020(a);
  803.     }
  804. }
  805.  
  806.  
  807. void bitrv2(int n, int *ip, double *a)
  808. {
  809.     int j, j1, k, k1, l, m, m2;
  810.     double xr, xi, yr, yi;
  811.     
  812.     ip[0] = 0;
  813.     l = n;
  814.     m = 1;
  815.     while ((m << 3) < l) {
  816.         l >>= 1;
  817.         for (j = 0; j < m; j++) {
  818.             ip[m + j] = ip[j] + l;
  819.         }
  820.         m <<= 1;
  821.     }
  822.     m2 = 2 * m;
  823.     if ((m << 3) == l) {
  824.         for (k = 0; k < m; k++) {
  825.             for (j = 0; j < k; j++) {
  826.                 j1 = 2 * j + ip[k];
  827.                 k1 = 2 * k + ip[j];
  828.                 xr = a[j1];
  829.                 xi = a[j1 + 1];
  830.                 yr = a[k1];
  831.                 yi = a[k1 + 1];
  832.                 a[j1] = yr;
  833.                 a[j1 + 1] = yi;
  834.                 a[k1] = xr;
  835.                 a[k1 + 1] = xi;
  836.                 j1 += m2;
  837.                 k1 += 2 * m2;
  838.                 xr = a[j1];
  839.                 xi = a[j1 + 1];
  840.                 yr = a[k1];
  841.                 yi = a[k1 + 1];
  842.                 a[j1] = yr;
  843.                 a[j1 + 1] = yi;
  844.                 a[k1] = xr;
  845.                 a[k1 + 1] = xi;
  846.                 j1 += m2;
  847.                 k1 -= m2;
  848.                 xr = a[j1];
  849.                 xi = a[j1 + 1];
  850.                 yr = a[k1];
  851.                 yi = a[k1 + 1];
  852.                 a[j1] = yr;
  853.                 a[j1 + 1] = yi;
  854.                 a[k1] = xr;
  855.                 a[k1 + 1] = xi;
  856.                 j1 += m2;
  857.                 k1 += 2 * m2;
  858.                 xr = a[j1];
  859.                 xi = a[j1 + 1];
  860.                 yr = a[k1];
  861.                 yi = a[k1 + 1];
  862.                 a[j1] = yr;
  863.                 a[j1 + 1] = yi;
  864.                 a[k1] = xr;
  865.                 a[k1 + 1] = xi;
  866.             }
  867.             j1 = 2 * k + m2 + ip[k];
  868.             k1 = j1 + m2;
  869.             xr = a[j1];
  870.             xi = a[j1 + 1];
  871.             yr = a[k1];
  872.             yi = a[k1 + 1];
  873.             a[j1] = yr;
  874.             a[j1 + 1] = yi;
  875.             a[k1] = xr;
  876.             a[k1 + 1] = xi;
  877.         }
  878.     } else {
  879.         for (k = 1; k < m; k++) {
  880.             for (j = 0; j < k; j++) {
  881.                 j1 = 2 * j + ip[k];
  882.                 k1 = 2 * k + ip[j];
  883.                 xr = a[j1];
  884.                 xi = a[j1 + 1];
  885.                 yr = a[k1];
  886.                 yi = a[k1 + 1];
  887.                 a[j1] = yr;
  888.                 a[j1 + 1] = yi;
  889.                 a[k1] = xr;
  890.                 a[k1 + 1] = xi;
  891.                 j1 += m2;
  892.                 k1 += m2;
  893.                 xr = a[j1];
  894.                 xi = a[j1 + 1];
  895.                 yr = a[k1];
  896.                 yi = a[k1 + 1];
  897.                 a[j1] = yr;
  898.                 a[j1 + 1] = yi;
  899.                 a[k1] = xr;
  900.                 a[k1 + 1] = xi;
  901.             }
  902.         }
  903.     }
  904. }
  905.  
  906.  
  907. void bitrv2conj(int n, int *ip, double *a)
  908. {
  909.     int j, j1, k, k1, l, m, m2;
  910.     double xr, xi, yr, yi;
  911.     
  912.     ip[0] = 0;
  913.     l = n;
  914.     m = 1;
  915.     while ((m << 3) < l) {
  916.         l >>= 1;
  917.         for (j = 0; j < m; j++) {
  918.             ip[m + j] = ip[j] + l;
  919.         }
  920.         m <<= 1;
  921.     }
  922.     m2 = 2 * m;
  923.     if ((m << 3) == l) {
  924.         for (k = 0; k < m; k++) {
  925.             for (j = 0; j < k; j++) {
  926.                 j1 = 2 * j + ip[k];
  927.                 k1 = 2 * k + ip[j];
  928.                 xr = a[j1];
  929.                 xi = -a[j1 + 1];
  930.                 yr = a[k1];
  931.                 yi = -a[k1 + 1];
  932.                 a[j1] = yr;
  933.                 a[j1 + 1] = yi;
  934.                 a[k1] = xr;
  935.                 a[k1 + 1] = xi;
  936.                 j1 += m2;
  937.                 k1 += 2 * m2;
  938.                 xr = a[j1];
  939.                 xi = -a[j1 + 1];
  940.                 yr = a[k1];
  941.                 yi = -a[k1 + 1];
  942.                 a[j1] = yr;
  943.                 a[j1 + 1] = yi;
  944.                 a[k1] = xr;
  945.                 a[k1 + 1] = xi;
  946.                 j1 += m2;
  947.                 k1 -= m2;
  948.                 xr = a[j1];
  949.                 xi = -a[j1 + 1];
  950.                 yr = a[k1];
  951.                 yi = -a[k1 + 1];
  952.                 a[j1] = yr;
  953.                 a[j1 + 1] = yi;
  954.                 a[k1] = xr;
  955.                 a[k1 + 1] = xi;
  956.                 j1 += m2;
  957.                 k1 += 2 * m2;
  958.                 xr = a[j1];
  959.                 xi = -a[j1 + 1];
  960.                 yr = a[k1];
  961.                 yi = -a[k1 + 1];
  962.                 a[j1] = yr;
  963.                 a[j1 + 1] = yi;
  964.                 a[k1] = xr;
  965.                 a[k1 + 1] = xi;
  966.             }
  967.             k1 = 2 * k + ip[k];
  968.             a[k1 + 1] = -a[k1 + 1];
  969.             j1 = k1 + m2;
  970.             k1 = j1 + m2;
  971.             xr = a[j1];
  972.             xi = -a[j1 + 1];
  973.             yr = a[k1];
  974.             yi = -a[k1 + 1];
  975.             a[j1] = yr;
  976.             a[j1 + 1] = yi;
  977.             a[k1] = xr;
  978.             a[k1 + 1] = xi;
  979.             k1 += m2;
  980.             a[k1 + 1] = -a[k1 + 1];
  981.         }
  982.     } else {
  983.         a[1] = -a[1];
  984.         a[m2 + 1] = -a[m2 + 1];
  985.         for (k = 1; k < m; k++) {
  986.             for (j = 0; j < k; j++) {
  987.                 j1 = 2 * j + ip[k];
  988.                 k1 = 2 * k + ip[j];
  989.                 xr = a[j1];
  990.                 xi = -a[j1 + 1];
  991.                 yr = a[k1];
  992.                 yi = -a[k1 + 1];
  993.                 a[j1] = yr;
  994.                 a[j1 + 1] = yi;
  995.                 a[k1] = xr;
  996.                 a[k1 + 1] = xi;
  997.                 j1 += m2;
  998.                 k1 += m2;
  999.                 xr = a[j1];
  1000.                 xi = -a[j1 + 1];
  1001.                 yr = a[k1];
  1002.                 yi = -a[k1 + 1];
  1003.                 a[j1] = yr;
  1004.                 a[j1 + 1] = yi;
  1005.                 a[k1] = xr;
  1006.                 a[k1 + 1] = xi;
  1007.             }
  1008.             k1 = 2 * k + ip[k];
  1009.             a[k1 + 1] = -a[k1 + 1];
  1010.             a[k1 + m2 + 1] = -a[k1 + m2 + 1];
  1011.         }
  1012.     }
  1013. }
  1014.  
  1015.  
  1016. void bitrv216(double *a)
  1017. {
  1018.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1019.         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
  1020.         x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
  1021.     
  1022.     x1r = a[2];
  1023.     x1i = a[3];
  1024.     x2r = a[4];
  1025.     x2i = a[5];
  1026.     x3r = a[6];
  1027.     x3i = a[7];
  1028.     x4r = a[8];
  1029.     x4i = a[9];
  1030.     x5r = a[10];
  1031.     x5i = a[11];
  1032.     x7r = a[14];
  1033.     x7i = a[15];
  1034.     x8r = a[16];
  1035.     x8i = a[17];
  1036.     x10r = a[20];
  1037.     x10i = a[21];
  1038.     x11r = a[22];
  1039.     x11i = a[23];
  1040.     x12r = a[24];
  1041.     x12i = a[25];
  1042.     x13r = a[26];
  1043.     x13i = a[27];
  1044.     x14r = a[28];
  1045.     x14i = a[29];
  1046.     a[2] = x8r;
  1047.     a[3] = x8i;
  1048.     a[4] = x4r;
  1049.     a[5] = x4i;
  1050.     a[6] = x12r;
  1051.     a[7] = x12i;
  1052.     a[8] = x2r;
  1053.     a[9] = x2i;
  1054.     a[10] = x10r;
  1055.     a[11] = x10i;
  1056.     a[14] = x14r;
  1057.     a[15] = x14i;
  1058.     a[16] = x1r;
  1059.     a[17] = x1i;
  1060.     a[20] = x5r;
  1061.     a[21] = x5i;
  1062.     a[22] = x13r;
  1063.     a[23] = x13i;
  1064.     a[24] = x3r;
  1065.     a[25] = x3i;
  1066.     a[26] = x11r;
  1067.     a[27] = x11i;
  1068.     a[28] = x7r;
  1069.     a[29] = x7i;
  1070. }
  1071.  
  1072.  
  1073. void bitrv216neg(double *a)
  1074. {
  1075.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1076.         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
  1077.         x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, 
  1078.         x13r, x13i, x14r, x14i, x15r, x15i;
  1079.     
  1080.     x1r = a[2];
  1081.     x1i = a[3];
  1082.     x2r = a[4];
  1083.     x2i = a[5];
  1084.     x3r = a[6];
  1085.     x3i = a[7];
  1086.     x4r = a[8];
  1087.     x4i = a[9];
  1088.     x5r = a[10];
  1089.     x5i = a[11];
  1090.     x6r = a[12];
  1091.     x6i = a[13];
  1092.     x7r = a[14];
  1093.     x7i = a[15];
  1094.     x8r = a[16];
  1095.     x8i = a[17];
  1096.     x9r = a[18];
  1097.     x9i = a[19];
  1098.     x10r = a[20];
  1099.     x10i = a[21];
  1100.     x11r = a[22];
  1101.     x11i = a[23];
  1102.     x12r = a[24];
  1103.     x12i = a[25];
  1104.     x13r = a[26];
  1105.     x13i = a[27];
  1106.     x14r = a[28];
  1107.     x14i = a[29];
  1108.     x15r = a[30];
  1109.     x15i = a[31];
  1110.     a[2] = x15r;
  1111.     a[3] = x15i;
  1112.     a[4] = x7r;
  1113.     a[5] = x7i;
  1114.     a[6] = x11r;
  1115.     a[7] = x11i;
  1116.     a[8] = x3r;
  1117.     a[9] = x3i;
  1118.     a[10] = x13r;
  1119.     a[11] = x13i;
  1120.     a[12] = x5r;
  1121.     a[13] = x5i;
  1122.     a[14] = x9r;
  1123.     a[15] = x9i;
  1124.     a[16] = x1r;
  1125.     a[17] = x1i;
  1126.     a[18] = x14r;
  1127.     a[19] = x14i;
  1128.     a[20] = x6r;
  1129.     a[21] = x6i;
  1130.     a[22] = x10r;
  1131.     a[23] = x10i;
  1132.     a[24] = x2r;
  1133.     a[25] = x2i;
  1134.     a[26] = x12r;
  1135.     a[27] = x12i;
  1136.     a[28] = x4r;
  1137.     a[29] = x4i;
  1138.     a[30] = x8r;
  1139.     a[31] = x8i;
  1140. }
  1141.  
  1142.  
  1143. void bitrv208(double *a)
  1144. {
  1145.     double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
  1146.     
  1147.     x1r = a[2];
  1148.     x1i = a[3];
  1149.     x3r = a[6];
  1150.     x3i = a[7];
  1151.     x4r = a[8];
  1152.     x4i = a[9];
  1153.     x6r = a[12];
  1154.     x6i = a[13];
  1155.     a[2] = x4r;
  1156.     a[3] = x4i;
  1157.     a[6] = x6r;
  1158.     a[7] = x6i;
  1159.     a[8] = x1r;
  1160.     a[9] = x1i;
  1161.     a[12] = x3r;
  1162.     a[13] = x3i;
  1163. }
  1164.  
  1165.  
  1166. void bitrv208neg(double *a)
  1167. {
  1168.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1169.         x5r, x5i, x6r, x6i, x7r, x7i;
  1170.     
  1171.     x1r = a[2];
  1172.     x1i = a[3];
  1173.     x2r = a[4];
  1174.     x2i = a[5];
  1175.     x3r = a[6];
  1176.     x3i = a[7];
  1177.     x4r = a[8];
  1178.     x4i = a[9];
  1179.     x5r = a[10];
  1180.     x5i = a[11];
  1181.     x6r = a[12];
  1182.     x6i = a[13];
  1183.     x7r = a[14];
  1184.     x7i = a[15];
  1185.     a[2] = x7r;
  1186.     a[3] = x7i;
  1187.     a[4] = x3r;
  1188.     a[5] = x3i;
  1189.     a[6] = x5r;
  1190.     a[7] = x5i;
  1191.     a[8] = x1r;
  1192.     a[9] = x1i;
  1193.     a[10] = x6r;
  1194.     a[11] = x6i;
  1195.     a[12] = x2r;
  1196.     a[13] = x2i;
  1197.     a[14] = x4r;
  1198.     a[15] = x4i;
  1199. }
  1200.  
  1201.  
  1202. void cftf1st(int n, double *a, double *w)
  1203. {
  1204.     int j, j0, j1, j2, j3, k, m, mh;
  1205.     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
  1206.         wd1r, wd1i, wd3r, wd3i;
  1207.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  1208.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
  1209.     
  1210.     mh = n >> 3;
  1211.     m = 2 * mh;
  1212.     j1 = m;
  1213.     j2 = j1 + m;
  1214.     j3 = j2 + m;
  1215.     x0r = a[0] + a[j2];
  1216.     x0i = a[1] + a[j2 + 1];
  1217.     x1r = a[0] - a[j2];
  1218.     x1i = a[1] - a[j2 + 1];
  1219.     x2r = a[j1] + a[j3];
  1220.     x2i = a[j1 + 1] + a[j3 + 1];
  1221.     x3r = a[j1] - a[j3];
  1222.     x3i = a[j1 + 1] - a[j3 + 1];
  1223.     a[0] = x0r + x2r;
  1224.     a[1] = x0i + x2i;
  1225.     a[j1] = x0r - x2r;
  1226.     a[j1 + 1] = x0i - x2i;
  1227.     a[j2] = x1r - x3i;
  1228.     a[j2 + 1] = x1i + x3r;
  1229.     a[j3] = x1r + x3i;
  1230.     a[j3 + 1] = x1i - x3r;
  1231.     wn4r = w[1];
  1232.     csc1 = w[2];
  1233.     csc3 = w[3];
  1234.     wd1r = 1;
  1235.     wd1i = 0;
  1236.     wd3r = 1;
  1237.     wd3i = 0;
  1238.     k = 0;
  1239.     for (j = 2; j < mh - 2; j += 4) {
  1240.         k += 4;
  1241.         wk1r = csc1 * (wd1r + w[k]);
  1242.         wk1i = csc1 * (wd1i + w[k + 1]);
  1243.         wk3r = csc3 * (wd3r + w[k + 2]);
  1244.         wk3i = csc3 * (wd3i - w[k + 3]);
  1245.         wd1r = w[k];
  1246.         wd1i = w[k + 1];
  1247.         wd3r = w[k + 2];
  1248.         wd3i = -w[k + 3];
  1249.         j1 = j + m;
  1250.         j2 = j1 + m;
  1251.         j3 = j2 + m;
  1252.         x0r = a[j] + a[j2];
  1253.         x0i = a[j + 1] + a[j2 + 1];
  1254.         x1r = a[j] - a[j2];
  1255.         x1i = a[j + 1] - a[j2 + 1];
  1256.         y0r = a[j + 2] + a[j2 + 2];
  1257.         y0i = a[j + 3] + a[j2 + 3];
  1258.         y1r = a[j + 2] - a[j2 + 2];
  1259.         y1i = a[j + 3] - a[j2 + 3];
  1260.         x2r = a[j1] + a[j3];
  1261.         x2i = a[j1 + 1] + a[j3 + 1];
  1262.         x3r = a[j1] - a[j3];
  1263.         x3i = a[j1 + 1] - a[j3 + 1];
  1264.         y2r = a[j1 + 2] + a[j3 + 2];
  1265.         y2i = a[j1 + 3] + a[j3 + 3];
  1266.         y3r = a[j1 + 2] - a[j3 + 2];
  1267.         y3i = a[j1 + 3] - a[j3 + 3];
  1268.         a[j] = x0r + x2r;
  1269.         a[j + 1] = x0i + x2i;
  1270.         a[j + 2] = y0r + y2r;
  1271.         a[j + 3] = y0i + y2i;
  1272.         a[j1] = x0r - x2r;
  1273.         a[j1 + 1] = x0i - x2i;
  1274.         a[j1 + 2] = y0r - y2r;
  1275.         a[j1 + 3] = y0i - y2i;
  1276.         x0r = x1r - x3i;
  1277.         x0i = x1i + x3r;
  1278.         a[j2] = wk1r * x0r - wk1i * x0i;
  1279.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  1280.         x0r = y1r - y3i;
  1281.         x0i = y1i + y3r;
  1282.         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
  1283.         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
  1284.         x0r = x1r + x3i;
  1285.         x0i = x1i - x3r;
  1286.         a[j3] = wk3r * x0r + wk3i * x0i;
  1287.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  1288.         x0r = y1r + y3i;
  1289.         x0i = y1i - y3r;
  1290.         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
  1291.         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
  1292.         j0 = m - j;
  1293.         j1 = j0 + m;
  1294.         j2 = j1 + m;
  1295.         j3 = j2 + m;
  1296.         x0r = a[j0] + a[j2];
  1297.         x0i = a[j0 + 1] + a[j2 + 1];
  1298.         x1r = a[j0] - a[j2];
  1299.         x1i = a[j0 + 1] - a[j2 + 1];
  1300.         y0r = a[j0 - 2] + a[j2 - 2];
  1301.         y0i = a[j0 - 1] + a[j2 - 1];
  1302.         y1r = a[j0 - 2] - a[j2 - 2];
  1303.         y1i = a[j0 - 1] - a[j2 - 1];
  1304.         x2r = a[j1] + a[j3];
  1305.         x2i = a[j1 + 1] + a[j3 + 1];
  1306.         x3r = a[j1] - a[j3];
  1307.         x3i = a[j1 + 1] - a[j3 + 1];
  1308.         y2r = a[j1 - 2] + a[j3 - 2];
  1309.         y2i = a[j1 - 1] + a[j3 - 1];
  1310.         y3r = a[j1 - 2] - a[j3 - 2];
  1311.         y3i = a[j1 - 1] - a[j3 - 1];
  1312.         a[j0] = x0r + x2r;
  1313.         a[j0 + 1] = x0i + x2i;
  1314.         a[j0 - 2] = y0r + y2r;
  1315.         a[j0 - 1] = y0i + y2i;
  1316.         a[j1] = x0r - x2r;
  1317.         a[j1 + 1] = x0i - x2i;
  1318.         a[j1 - 2] = y0r - y2r;
  1319.         a[j1 - 1] = y0i - y2i;
  1320.         x0r = x1r - x3i;
  1321.         x0i = x1i + x3r;
  1322.         a[j2] = wk1i * x0r - wk1r * x0i;
  1323.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  1324.         x0r = y1r - y3i;
  1325.         x0i = y1i + y3r;
  1326.         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
  1327.         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
  1328.         x0r = x1r + x3i;
  1329.         x0i = x1i - x3r;
  1330.         a[j3] = wk3i * x0r + wk3r * x0i;
  1331.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  1332.         x0r = y1r + y3i;
  1333.         x0i = y1i - y3r;
  1334.         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
  1335.         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
  1336.     }
  1337.     wk1r = csc1 * (wd1r + wn4r);
  1338.     wk1i = csc1 * (wd1i + wn4r);
  1339.     wk3r = csc3 * (wd3r - wn4r);
  1340.     wk3i = csc3 * (wd3i - wn4r);
  1341.     j0 = mh;
  1342.     j1 = j0 + m;
  1343.     j2 = j1 + m;
  1344.     j3 = j2 + m;
  1345.     x0r = a[j0 - 2] + a[j2 - 2];
  1346.     x0i = a[j0 - 1] + a[j2 - 1];
  1347.     x1r = a[j0 - 2] - a[j2 - 2];
  1348.     x1i = a[j0 - 1] - a[j2 - 1];
  1349.     x2r = a[j1 - 2] + a[j3 - 2];
  1350.     x2i = a[j1 - 1] + a[j3 - 1];
  1351.     x3r = a[j1 - 2] - a[j3 - 2];
  1352.     x3i = a[j1 - 1] - a[j3 - 1];
  1353.     a[j0 - 2] = x0r + x2r;
  1354.     a[j0 - 1] = x0i + x2i;
  1355.     a[j1 - 2] = x0r - x2r;
  1356.     a[j1 - 1] = x0i - x2i;
  1357.     x0r = x1r - x3i;
  1358.     x0i = x1i + x3r;
  1359.     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
  1360.     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
  1361.     x0r = x1r + x3i;
  1362.     x0i = x1i - x3r;
  1363.     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
  1364.     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
  1365.     x0r = a[j0] + a[j2];
  1366.     x0i = a[j0 + 1] + a[j2 + 1];
  1367.     x1r = a[j0] - a[j2];
  1368.     x1i = a[j0 + 1] - a[j2 + 1];
  1369.     x2r = a[j1] + a[j3];
  1370.     x2i = a[j1 + 1] + a[j3 + 1];
  1371.     x3r = a[j1] - a[j3];
  1372.     x3i = a[j1 + 1] - a[j3 + 1];
  1373.     a[j0] = x0r + x2r;
  1374.     a[j0 + 1] = x0i + x2i;
  1375.     a[j1] = x0r - x2r;
  1376.     a[j1 + 1] = x0i - x2i;
  1377.     x0r = x1r - x3i;
  1378.     x0i = x1i + x3r;
  1379.     a[j2] = wn4r * (x0r - x0i);
  1380.     a[j2 + 1] = wn4r * (x0i + x0r);
  1381.     x0r = x1r + x3i;
  1382.     x0i = x1i - x3r;
  1383.     a[j3] = -wn4r * (x0r + x0i);
  1384.     a[j3 + 1] = -wn4r * (x0i - x0r);
  1385.     x0r = a[j0 + 2] + a[j2 + 2];
  1386.     x0i = a[j0 + 3] + a[j2 + 3];
  1387.     x1r = a[j0 + 2] - a[j2 + 2];
  1388.     x1i = a[j0 + 3] - a[j2 + 3];
  1389.     x2r = a[j1 + 2] + a[j3 + 2];
  1390.     x2i = a[j1 + 3] + a[j3 + 3];
  1391.     x3r = a[j1 + 2] - a[j3 + 2];
  1392.     x3i = a[j1 + 3] - a[j3 + 3];
  1393.     a[j0 + 2] = x0r + x2r;
  1394.     a[j0 + 3] = x0i + x2i;
  1395.     a[j1 + 2] = x0r - x2r;
  1396.     a[j1 + 3] = x0i - x2i;
  1397.     x0r = x1r - x3i;
  1398.     x0i = x1i + x3r;
  1399.     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
  1400.     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
  1401.     x0r = x1r + x3i;
  1402.     x0i = x1i - x3r;
  1403.     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
  1404.     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
  1405. }
  1406.  
  1407.  
  1408. void cftb1st(int n, double *a, double *w)
  1409. {
  1410.     int j, j0, j1, j2, j3, k, m, mh;
  1411.     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
  1412.         wd1r, wd1i, wd3r, wd3i;
  1413.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  1414.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
  1415.     
  1416.     mh = n >> 3;
  1417.     m = 2 * mh;
  1418.     j1 = m;
  1419.     j2 = j1 + m;
  1420.     j3 = j2 + m;
  1421.     x0r = a[0] + a[j2];
  1422.     x0i = -a[1] - a[j2 + 1];
  1423.     x1r = a[0] - a[j2];
  1424.     x1i = -a[1] + a[j2 + 1];
  1425.     x2r = a[j1] + a[j3];
  1426.     x2i = a[j1 + 1] + a[j3 + 1];
  1427.     x3r = a[j1] - a[j3];
  1428.     x3i = a[j1 + 1] - a[j3 + 1];
  1429.     a[0] = x0r + x2r;
  1430.     a[1] = x0i - x2i;
  1431.     a[j1] = x0r - x2r;
  1432.     a[j1 + 1] = x0i + x2i;
  1433.     a[j2] = x1r + x3i;
  1434.     a[j2 + 1] = x1i + x3r;
  1435.     a[j3] = x1r - x3i;
  1436.     a[j3 + 1] = x1i - x3r;
  1437.     wn4r = w[1];
  1438.     csc1 = w[2];
  1439.     csc3 = w[3];
  1440.     wd1r = 1;
  1441.     wd1i = 0;
  1442.     wd3r = 1;
  1443.     wd3i = 0;
  1444.     k = 0;
  1445.     for (j = 2; j < mh - 2; j += 4) {
  1446.         k += 4;
  1447.         wk1r = csc1 * (wd1r + w[k]);
  1448.         wk1i = csc1 * (wd1i + w[k + 1]);
  1449.         wk3r = csc3 * (wd3r + w[k + 2]);
  1450.         wk3i = csc3 * (wd3i - w[k + 3]);
  1451.         wd1r = w[k];
  1452.         wd1i = w[k + 1];
  1453.         wd3r = w[k + 2];
  1454.         wd3i = -w[k + 3];
  1455.         j1 = j + m;
  1456.         j2 = j1 + m;
  1457.         j3 = j2 + m;
  1458.         x0r = a[j] + a[j2];
  1459.         x0i = -a[j + 1] - a[j2 + 1];
  1460.         x1r = a[j] - a[j2];
  1461.         x1i = -a[j + 1] + a[j2 + 1];
  1462.         y0r = a[j + 2] + a[j2 + 2];
  1463.         y0i = -a[j + 3] - a[j2 + 3];
  1464.         y1r = a[j + 2] - a[j2 + 2];
  1465.         y1i = -a[j + 3] + a[j2 + 3];
  1466.         x2r = a[j1] + a[j3];
  1467.         x2i = a[j1 + 1] + a[j3 + 1];
  1468.         x3r = a[j1] - a[j3];
  1469.         x3i = a[j1 + 1] - a[j3 + 1];
  1470.         y2r = a[j1 + 2] + a[j3 + 2];
  1471.         y2i = a[j1 + 3] + a[j3 + 3];
  1472.         y3r = a[j1 + 2] - a[j3 + 2];
  1473.         y3i = a[j1 + 3] - a[j3 + 3];
  1474.         a[j] = x0r + x2r;
  1475.         a[j + 1] = x0i - x2i;
  1476.         a[j + 2] = y0r + y2r;
  1477.         a[j + 3] = y0i - y2i;
  1478.         a[j1] = x0r - x2r;
  1479.         a[j1 + 1] = x0i + x2i;
  1480.         a[j1 + 2] = y0r - y2r;
  1481.         a[j1 + 3] = y0i + y2i;
  1482.         x0r = x1r + x3i;
  1483.         x0i = x1i + x3r;
  1484.         a[j2] = wk1r * x0r - wk1i * x0i;
  1485.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  1486.         x0r = y1r + y3i;
  1487.         x0i = y1i + y3r;
  1488.         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
  1489.         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
  1490.         x0r = x1r - x3i;
  1491.         x0i = x1i - x3r;
  1492.         a[j3] = wk3r * x0r + wk3i * x0i;
  1493.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  1494.         x0r = y1r - y3i;
  1495.         x0i = y1i - y3r;
  1496.         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
  1497.         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
  1498.         j0 = m - j;
  1499.         j1 = j0 + m;
  1500.         j2 = j1 + m;
  1501.         j3 = j2 + m;
  1502.         x0r = a[j0] + a[j2];
  1503.         x0i = -a[j0 + 1] - a[j2 + 1];
  1504.         x1r = a[j0] - a[j2];
  1505.         x1i = -a[j0 + 1] + a[j2 + 1];
  1506.         y0r = a[j0 - 2] + a[j2 - 2];
  1507.         y0i = -a[j0 - 1] - a[j2 - 1];
  1508.         y1r = a[j0 - 2] - a[j2 - 2];
  1509.         y1i = -a[j0 - 1] + a[j2 - 1];
  1510.         x2r = a[j1] + a[j3];
  1511.         x2i = a[j1 + 1] + a[j3 + 1];
  1512.         x3r = a[j1] - a[j3];
  1513.         x3i = a[j1 + 1] - a[j3 + 1];
  1514.         y2r = a[j1 - 2] + a[j3 - 2];
  1515.         y2i = a[j1 - 1] + a[j3 - 1];
  1516.         y3r = a[j1 - 2] - a[j3 - 2];
  1517.         y3i = a[j1 - 1] - a[j3 - 1];
  1518.         a[j0] = x0r + x2r;
  1519.         a[j0 + 1] = x0i - x2i;
  1520.         a[j0 - 2] = y0r + y2r;
  1521.         a[j0 - 1] = y0i - y2i;
  1522.         a[j1] = x0r - x2r;
  1523.         a[j1 + 1] = x0i + x2i;
  1524.         a[j1 - 2] = y0r - y2r;
  1525.         a[j1 - 1] = y0i + y2i;
  1526.         x0r = x1r + x3i;
  1527.         x0i = x1i + x3r;
  1528.         a[j2] = wk1i * x0r - wk1r * x0i;
  1529.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  1530.         x0r = y1r + y3i;
  1531.         x0i = y1i + y3r;
  1532.         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
  1533.         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
  1534.         x0r = x1r - x3i;
  1535.         x0i = x1i - x3r;
  1536.         a[j3] = wk3i * x0r + wk3r * x0i;
  1537.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  1538.         x0r = y1r - y3i;
  1539.         x0i = y1i - y3r;
  1540.         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
  1541.         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
  1542.     }
  1543.     wk1r = csc1 * (wd1r + wn4r);
  1544.     wk1i = csc1 * (wd1i + wn4r);
  1545.     wk3r = csc3 * (wd3r - wn4r);
  1546.     wk3i = csc3 * (wd3i - wn4r);
  1547.     j0 = mh;
  1548.     j1 = j0 + m;
  1549.     j2 = j1 + m;
  1550.     j3 = j2 + m;
  1551.     x0r = a[j0 - 2] + a[j2 - 2];
  1552.     x0i = -a[j0 - 1] - a[j2 - 1];
  1553.     x1r = a[j0 - 2] - a[j2 - 2];
  1554.     x1i = -a[j0 - 1] + a[j2 - 1];
  1555.     x2r = a[j1 - 2] + a[j3 - 2];
  1556.     x2i = a[j1 - 1] + a[j3 - 1];
  1557.     x3r = a[j1 - 2] - a[j3 - 2];
  1558.     x3i = a[j1 - 1] - a[j3 - 1];
  1559.     a[j0 - 2] = x0r + x2r;
  1560.     a[j0 - 1] = x0i - x2i;
  1561.     a[j1 - 2] = x0r - x2r;
  1562.     a[j1 - 1] = x0i + x2i;
  1563.     x0r = x1r + x3i;
  1564.     x0i = x1i + x3r;
  1565.     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
  1566.     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
  1567.     x0r = x1r - x3i;
  1568.     x0i = x1i - x3r;
  1569.     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
  1570.     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
  1571.     x0r = a[j0] + a[j2];
  1572.     x0i = -a[j0 + 1] - a[j2 + 1];
  1573.     x1r = a[j0] - a[j2];
  1574.     x1i = -a[j0 + 1] + a[j2 + 1];
  1575.     x2r = a[j1] + a[j3];
  1576.     x2i = a[j1 + 1] + a[j3 + 1];
  1577.     x3r = a[j1] - a[j3];
  1578.     x3i = a[j1 + 1] - a[j3 + 1];
  1579.     a[j0] = x0r + x2r;
  1580.     a[j0 + 1] = x0i - x2i;
  1581.     a[j1] = x0r - x2r;
  1582.     a[j1 + 1] = x0i + x2i;
  1583.     x0r = x1r + x3i;
  1584.     x0i = x1i + x3r;
  1585.     a[j2] = wn4r * (x0r - x0i);
  1586.     a[j2 + 1] = wn4r * (x0i + x0r);
  1587.     x0r = x1r - x3i;
  1588.     x0i = x1i - x3r;
  1589.     a[j3] = -wn4r * (x0r + x0i);
  1590.     a[j3 + 1] = -wn4r * (x0i - x0r);
  1591.     x0r = a[j0 + 2] + a[j2 + 2];
  1592.     x0i = -a[j0 + 3] - a[j2 + 3];
  1593.     x1r = a[j0 + 2] - a[j2 + 2];
  1594.     x1i = -a[j0 + 3] + a[j2 + 3];
  1595.     x2r = a[j1 + 2] + a[j3 + 2];
  1596.     x2i = a[j1 + 3] + a[j3 + 3];
  1597.     x3r = a[j1 + 2] - a[j3 + 2];
  1598.     x3i = a[j1 + 3] - a[j3 + 3];
  1599.     a[j0 + 2] = x0r + x2r;
  1600.     a[j0 + 3] = x0i - x2i;
  1601.     a[j1 + 2] = x0r - x2r;
  1602.     a[j1 + 3] = x0i + x2i;
  1603.     x0r = x1r + x3i;
  1604.     x0i = x1i + x3r;
  1605.     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
  1606.     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
  1607.     x0r = x1r - x3i;
  1608.     x0i = x1i - x3r;
  1609.     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
  1610.     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
  1611. }
  1612.  
  1613.  
  1614. void cftrec1(int n, double *a, int nw, double *w)
  1615. {
  1616.     void cftrec1(int n, double *a, int nw, double *w);
  1617.     void cftrec2(int n, double *a, int nw, double *w);
  1618.     void cftmdl1(int n, double *a, double *w);
  1619.     void cftexp1(int n, double *a, int nw, double *w);
  1620.     int m;
  1621.     
  1622.     m = n >> 2;
  1623.     cftmdl1(n, a, &w[nw - 2 * m]);
  1624.     if (n > CDFT_RECURSIVE_N) {
  1625.         cftrec1(m, a, nw, w);
  1626.         cftrec2(m, &a[m], nw, w);
  1627.         cftrec1(m, &a[2 * m], nw, w);
  1628.         cftrec1(m, &a[3 * m], nw, w);
  1629.     } else {
  1630.         cftexp1(n, a, nw, w);
  1631.     }
  1632. }
  1633.  
  1634.  
  1635. void cftrec2(int n, double *a, int nw, double *w)
  1636. {
  1637.     void cftrec1(int n, double *a, int nw, double *w);
  1638.     void cftrec2(int n, double *a, int nw, double *w);
  1639.     void cftmdl2(int n, double *a, double *w);
  1640.     void cftexp2(int n, double *a, int nw, double *w);
  1641.     int m;
  1642.     
  1643.     m = n >> 2;
  1644.     cftmdl2(n, a, &w[nw - n]);
  1645.     if (n > CDFT_RECURSIVE_N) {
  1646.         cftrec1(m, a, nw, w);
  1647.         cftrec2(m, &a[m], nw, w);
  1648.         cftrec1(m, &a[2 * m], nw, w);
  1649.         cftrec2(m, &a[3 * m], nw, w);
  1650.     } else {
  1651.         cftexp2(n, a, nw, w);
  1652.     }
  1653. }
  1654.  
  1655.  
  1656. void cftexp1(int n, double *a, int nw, double *w)
  1657. {
  1658.     void cftmdl1(int n, double *a, double *w);
  1659.     void cftmdl2(int n, double *a, double *w);
  1660.     void cftfx41(int n, double *a, int nw, double *w);
  1661.     void cftfx42(int n, double *a, int nw, double *w);
  1662.     int j, k, l;
  1663.     
  1664.     l = n >> 2;
  1665.     while (l > 128) {
  1666.         for (k = l; k < n; k <<= 2) {
  1667.             for (j = k - l; j < n; j += 4 * k) {
  1668.                 cftmdl1(l, &a[j], &w[nw - (l >> 1)]);
  1669.                 cftmdl2(l, &a[k + j], &w[nw - l]);
  1670.                 cftmdl1(l, &a[2 * k + j], &w[nw - (l >> 1)]);
  1671.             }
  1672.         }
  1673.         cftmdl1(l, &a[n - l], &w[nw - (l >> 1)]);
  1674.         l >>= 2;
  1675.     }
  1676.     for (k = l; k < n; k <<= 2) {
  1677.         for (j = k - l; j < n; j += 4 * k) {
  1678.             cftmdl1(l, &a[j], &w[nw - (l >> 1)]);
  1679.             cftfx41(l, &a[j], nw, w);
  1680.             cftmdl2(l, &a[k + j], &w[nw - l]);
  1681.             cftfx42(l, &a[k + j], nw, w);
  1682.             cftmdl1(l, &a[2 * k + j], &w[nw - (l >> 1)]);
  1683.             cftfx41(l, &a[2 * k + j], nw, w);
  1684.         }
  1685.     }
  1686.     cftmdl1(l, &a[n - l], &w[nw - (l >> 1)]);
  1687.     cftfx41(l, &a[n - l], nw, w);
  1688. }
  1689.  
  1690.  
  1691. void cftexp2(int n, double *a, int nw, double *w)
  1692. {
  1693.     void cftmdl1(int n, double *a, double *w);
  1694.     void cftmdl2(int n, double *a, double *w);
  1695.     void cftfx41(int n, double *a, int nw, double *w);
  1696.     void cftfx42(int n, double *a, int nw, double *w);
  1697.     int j, k, l, m;
  1698.     
  1699.     m = n >> 1;
  1700.     l = n >> 2;
  1701.     while (l > 128) {
  1702.         for (k = l; k < m; k <<= 2) {
  1703.             for (j = k - l; j < m; j += 2 * k) {
  1704.                 cftmdl1(l, &a[j], &w[nw - (l >> 1)]);
  1705.                 cftmdl1(l, &a[m + j], &w[nw - (l >> 1)]);
  1706.             }
  1707.             for (j = 2 * k - l; j < m; j += 4 * k) {
  1708.                 cftmdl2(l, &a[j], &w[nw - l]);
  1709.                 cftmdl2(l, &a[m + j], &w[nw - l]);
  1710.             }
  1711.         }
  1712.         l >>= 2;
  1713.     }
  1714.     for (k = l; k < m; k <<= 2) {
  1715.         for (j = k - l; j < m; j += 2 * k) {
  1716.             cftmdl1(l, &a[j], &w[nw - (l >> 1)]);
  1717.             cftfx41(l, &a[j], nw, w);
  1718.             cftmdl1(l, &a[m + j], &w[nw - (l >> 1)]);
  1719.             cftfx41(l, &a[m + j], nw, w);
  1720.         }
  1721.         for (j = 2 * k - l; j < m; j += 4 * k) {
  1722.             cftmdl2(l, &a[j], &w[nw - l]);
  1723.             cftfx42(l, &a[j], nw, w);
  1724.             cftmdl2(l, &a[m + j], &w[nw - l]);
  1725.             cftfx42(l, &a[m + j], nw, w);
  1726.         }
  1727.     }
  1728. }
  1729.  
  1730.  
  1731. void cftmdl1(int n, double *a, double *w)
  1732. {
  1733.     int j, j0, j1, j2, j3, k, m, mh;
  1734.     double wn4r, wk1r, wk1i, wk3r, wk3i;
  1735.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  1736.     
  1737.     mh = n >> 3;
  1738.     m = 2 * mh;
  1739.     j1 = m;
  1740.     j2 = j1 + m;
  1741.     j3 = j2 + m;
  1742.     x0r = a[0] + a[j2];
  1743.     x0i = a[1] + a[j2 + 1];
  1744.     x1r = a[0] - a[j2];
  1745.     x1i = a[1] - a[j2 + 1];
  1746.     x2r = a[j1] + a[j3];
  1747.     x2i = a[j1 + 1] + a[j3 + 1];
  1748.     x3r = a[j1] - a[j3];
  1749.     x3i = a[j1 + 1] - a[j3 + 1];
  1750.     a[0] = x0r + x2r;
  1751.     a[1] = x0i + x2i;
  1752.     a[j1] = x0r - x2r;
  1753.     a[j1 + 1] = x0i - x2i;
  1754.     a[j2] = x1r - x3i;
  1755.     a[j2 + 1] = x1i + x3r;
  1756.     a[j3] = x1r + x3i;
  1757.     a[j3 + 1] = x1i - x3r;
  1758.     wn4r = w[1];
  1759.     k = 0;
  1760.     for (j = 2; j < mh; j += 2) {
  1761.         k += 4;
  1762.         wk1r = w[k];
  1763.         wk1i = w[k + 1];
  1764.         wk3r = w[k + 2];
  1765.         wk3i = -w[k + 3];
  1766.         j1 = j + m;
  1767.         j2 = j1 + m;
  1768.         j3 = j2 + m;
  1769.         x0r = a[j] + a[j2];
  1770.         x0i = a[j + 1] + a[j2 + 1];
  1771.         x1r = a[j] - a[j2];
  1772.         x1i = a[j + 1] - a[j2 + 1];
  1773.         x2r = a[j1] + a[j3];
  1774.         x2i = a[j1 + 1] + a[j3 + 1];
  1775.         x3r = a[j1] - a[j3];
  1776.         x3i = a[j1 + 1] - a[j3 + 1];
  1777.         a[j] = x0r + x2r;
  1778.         a[j + 1] = x0i + x2i;
  1779.         a[j1] = x0r - x2r;
  1780.         a[j1 + 1] = x0i - x2i;
  1781.         x0r = x1r - x3i;
  1782.         x0i = x1i + x3r;
  1783.         a[j2] = wk1r * x0r - wk1i * x0i;
  1784.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  1785.         x0r = x1r + x3i;
  1786.         x0i = x1i - x3r;
  1787.         a[j3] = wk3r * x0r + wk3i * x0i;
  1788.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  1789.         j0 = m - j;
  1790.         j1 = j0 + m;
  1791.         j2 = j1 + m;
  1792.         j3 = j2 + m;
  1793.         x0r = a[j0] + a[j2];
  1794.         x0i = a[j0 + 1] + a[j2 + 1];
  1795.         x1r = a[j0] - a[j2];
  1796.         x1i = a[j0 + 1] - a[j2 + 1];
  1797.         x2r = a[j1] + a[j3];
  1798.         x2i = a[j1 + 1] + a[j3 + 1];
  1799.         x3r = a[j1] - a[j3];
  1800.         x3i = a[j1 + 1] - a[j3 + 1];
  1801.         a[j0] = x0r + x2r;
  1802.         a[j0 + 1] = x0i + x2i;
  1803.         a[j1] = x0r - x2r;
  1804.         a[j1 + 1] = x0i - x2i;
  1805.         x0r = x1r - x3i;
  1806.         x0i = x1i + x3r;
  1807.         a[j2] = wk1i * x0r - wk1r * x0i;
  1808.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  1809.         x0r = x1r + x3i;
  1810.         x0i = x1i - x3r;
  1811.         a[j3] = wk3i * x0r + wk3r * x0i;
  1812.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  1813.     }
  1814.     j0 = mh;
  1815.     j1 = j0 + m;
  1816.     j2 = j1 + m;
  1817.     j3 = j2 + m;
  1818.     x0r = a[j0] + a[j2];
  1819.     x0i = a[j0 + 1] + a[j2 + 1];
  1820.     x1r = a[j0] - a[j2];
  1821.     x1i = a[j0 + 1] - a[j2 + 1];
  1822.     x2r = a[j1] + a[j3];
  1823.     x2i = a[j1 + 1] + a[j3 + 1];
  1824.     x3r = a[j1] - a[j3];
  1825.     x3i = a[j1 + 1] - a[j3 + 1];
  1826.     a[j0] = x0r + x2r;
  1827.     a[j0 + 1] = x0i + x2i;
  1828.     a[j1] = x0r - x2r;
  1829.     a[j1 + 1] = x0i - x2i;
  1830.     x0r = x1r - x3i;
  1831.     x0i = x1i + x3r;
  1832.     a[j2] = wn4r * (x0r - x0i);
  1833.     a[j2 + 1] = wn4r * (x0i + x0r);
  1834.     x0r = x1r + x3i;
  1835.     x0i = x1i - x3r;
  1836.     a[j3] = -wn4r * (x0r + x0i);
  1837.     a[j3 + 1] = -wn4r * (x0i - x0r);
  1838. }
  1839.  
  1840.  
  1841. void cftmdl2(int n, double *a, double *w)
  1842. {
  1843.     int j, j0, j1, j2, j3, k, kr, m, mh;
  1844.     double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
  1845.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
  1846.     
  1847.     mh = n >> 3;
  1848.     m = 2 * mh;
  1849.     wn4r = w[1];
  1850.     j1 = m;
  1851.     j2 = j1 + m;
  1852.     j3 = j2 + m;
  1853.     x0r = a[0] - a[j2 + 1];
  1854.     x0i = a[1] + a[j2];
  1855.     x1r = a[0] + a[j2 + 1];
  1856.     x1i = a[1] - a[j2];
  1857.     x2r = a[j1] - a[j3 + 1];
  1858.     x2i = a[j1 + 1] + a[j3];
  1859.     x3r = a[j1] + a[j3 + 1];
  1860.     x3i = a[j1 + 1] - a[j3];
  1861.     y0r = wn4r * (x2r - x2i);
  1862.     y0i = wn4r * (x2i + x2r);
  1863.     a[0] = x0r + y0r;
  1864.     a[1] = x0i + y0i;
  1865.     a[j1] = x0r - y0r;
  1866.     a[j1 + 1] = x0i - y0i;
  1867.     y0r = wn4r * (x3r - x3i);
  1868.     y0i = wn4r * (x3i + x3r);
  1869.     a[j2] = x1r - y0i;
  1870.     a[j2 + 1] = x1i + y0r;
  1871.     a[j3] = x1r + y0i;
  1872.     a[j3 + 1] = x1i - y0r;
  1873.     k = 0;
  1874.     kr = 2 * m;
  1875.     for (j = 2; j < mh; j += 2) {
  1876.         k += 4;
  1877.         wk1r = w[k];
  1878.         wk1i = w[k + 1];
  1879.         wk3r = w[k + 2];
  1880.         wk3i = -w[k + 3];
  1881.         kr -= 4;
  1882.         wd1i = w[kr];
  1883.         wd1r = w[kr + 1];
  1884.         wd3i = w[kr + 2];
  1885.         wd3r = -w[kr + 3];
  1886.         j1 = j + m;
  1887.         j2 = j1 + m;
  1888.         j3 = j2 + m;
  1889.         x0r = a[j] - a[j2 + 1];
  1890.         x0i = a[j + 1] + a[j2];
  1891.         x1r = a[j] + a[j2 + 1];
  1892.         x1i = a[j + 1] - a[j2];
  1893.         x2r = a[j1] - a[j3 + 1];
  1894.         x2i = a[j1 + 1] + a[j3];
  1895.         x3r = a[j1] + a[j3 + 1];
  1896.         x3i = a[j1 + 1] - a[j3];
  1897.         y0r = wk1r * x0r - wk1i * x0i;
  1898.         y0i = wk1r * x0i + wk1i * x0r;
  1899.         y2r = wd1r * x2r - wd1i * x2i;
  1900.         y2i = wd1r * x2i + wd1i * x2r;
  1901.         a[j] = y0r + y2r;
  1902.         a[j + 1] = y0i + y2i;
  1903.         a[j1] = y0r - y2r;
  1904.         a[j1 + 1] = y0i - y2i;
  1905.         y0r = wk3r * x1r + wk3i * x1i;
  1906.         y0i = wk3r * x1i - wk3i * x1r;
  1907.         y2r = wd3r * x3r + wd3i * x3i;
  1908.         y2i = wd3r * x3i - wd3i * x3r;
  1909.         a[j2] = y0r + y2r;
  1910.         a[j2 + 1] = y0i + y2i;
  1911.         a[j3] = y0r - y2r;
  1912.         a[j3 + 1] = y0i - y2i;
  1913.         j0 = m - j;
  1914.         j1 = j0 + m;
  1915.         j2 = j1 + m;
  1916.         j3 = j2 + m;
  1917.         x0r = a[j0] - a[j2 + 1];
  1918.         x0i = a[j0 + 1] + a[j2];
  1919.         x1r = a[j0] + a[j2 + 1];
  1920.         x1i = a[j0 + 1] - a[j2];
  1921.         x2r = a[j1] - a[j3 + 1];
  1922.         x2i = a[j1 + 1] + a[j3];
  1923.         x3r = a[j1] + a[j3 + 1];
  1924.         x3i = a[j1 + 1] - a[j3];
  1925.         y0r = wd1i * x0r - wd1r * x0i;
  1926.         y0i = wd1i * x0i + wd1r * x0r;
  1927.         y2r = wk1i * x2r - wk1r * x2i;
  1928.         y2i = wk1i * x2i + wk1r * x2r;
  1929.         a[j0] = y0r + y2r;
  1930.         a[j0 + 1] = y0i + y2i;
  1931.         a[j1] = y0r - y2r;
  1932.         a[j1 + 1] = y0i - y2i;
  1933.         y0r = wd3i * x1r + wd3r * x1i;
  1934.         y0i = wd3i * x1i - wd3r * x1r;
  1935.         y2r = wk3i * x3r + wk3r * x3i;
  1936.         y2i = wk3i * x3i - wk3r * x3r;
  1937.         a[j2] = y0r + y2r;
  1938.         a[j2 + 1] = y0i + y2i;
  1939.         a[j3] = y0r - y2r;
  1940.         a[j3 + 1] = y0i - y2i;
  1941.     }
  1942.     wk1r = w[m];
  1943.     wk1i = w[m + 1];
  1944.     j0 = mh;
  1945.     j1 = j0 + m;
  1946.     j2 = j1 + m;
  1947.     j3 = j2 + m;
  1948.     x0r = a[j0] - a[j2 + 1];
  1949.     x0i = a[j0 + 1] + a[j2];
  1950.     x1r = a[j0] + a[j2 + 1];
  1951.     x1i = a[j0 + 1] - a[j2];
  1952.     x2r = a[j1] - a[j3 + 1];
  1953.     x2i = a[j1 + 1] + a[j3];
  1954.     x3r = a[j1] + a[j3 + 1];
  1955.     x3i = a[j1 + 1] - a[j3];
  1956.     y0r = wk1r * x0r - wk1i * x0i;
  1957.     y0i = wk1r * x0i + wk1i * x0r;
  1958.     y2r = wk1i * x2r - wk1r * x2i;
  1959.     y2i = wk1i * x2i + wk1r * x2r;
  1960.     a[j0] = y0r + y2r;
  1961.     a[j0 + 1] = y0i + y2i;
  1962.     a[j1] = y0r - y2r;
  1963.     a[j1 + 1] = y0i - y2i;
  1964.     y0r = wk1i * x1r - wk1r * x1i;
  1965.     y0i = wk1i * x1i + wk1r * x1r;
  1966.     y2r = wk1r * x3r - wk1i * x3i;
  1967.     y2i = wk1r * x3i + wk1i * x3r;
  1968.     a[j2] = y0r - y2r;
  1969.     a[j2 + 1] = y0i - y2i;
  1970.     a[j3] = y0r + y2r;
  1971.     a[j3 + 1] = y0i + y2i;
  1972. }
  1973.  
  1974.  
  1975. void cftfx41(int n, double *a, int nw, double *w)
  1976. {
  1977.     void cftf161(double *a, double *w);
  1978.     void cftf162(double *a, double *w);
  1979.     void cftf081(double *a, double *w);
  1980.     void cftf082(double *a, double *w);
  1981.     
  1982.     if (n == 128) {
  1983.         cftf161(a, &w[nw - 8]);
  1984.         cftf162(&a[32], &w[nw - 32]);
  1985.         cftf161(&a[64], &w[nw - 8]);
  1986.         cftf161(&a[96], &w[nw - 8]);
  1987.     } else {
  1988.         cftf081(a, &w[nw - 16]);
  1989.         cftf082(&a[16], &w[nw - 16]);
  1990.         cftf081(&a[32], &w[nw - 16]);
  1991.         cftf081(&a[48], &w[nw - 16]);
  1992.     }
  1993. }
  1994.  
  1995.  
  1996. void cftfx42(int n, double *a, int nw, double *w)
  1997. {
  1998.     void cftf161(double *a, double *w);
  1999.     void cftf162(double *a, double *w);
  2000.     void cftf081(double *a, double *w);
  2001.     void cftf082(double *a, double *w);
  2002.     
  2003.     if (n == 128) {
  2004.         cftf161(a, &w[nw - 8]);
  2005.         cftf162(&a[32], &w[nw - 32]);
  2006.         cftf161(&a[64], &w[nw - 8]);
  2007.         cftf162(&a[96], &w[nw - 32]);
  2008.     } else {
  2009.         cftf081(a, &w[nw - 16]);
  2010.         cftf082(&a[16], &w[nw - 16]);
  2011.         cftf081(&a[32], &w[nw - 16]);
  2012.         cftf082(&a[48], &w[nw - 16]);
  2013.     }
  2014. }
  2015.  
  2016.  
  2017. void cftf161(double *a, double *w)
  2018. {
  2019.     double wn4r, wk1r, wk1i, 
  2020.         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  2021.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2022.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
  2023.         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
  2024.         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
  2025.     
  2026.     wn4r = w[1];
  2027.     wk1i = wn4r * w[2];
  2028.     wk1r = wk1i + w[2];
  2029.     x0r = a[0] + a[16];
  2030.     x0i = a[1] + a[17];
  2031.     x1r = a[0] - a[16];
  2032.     x1i = a[1] - a[17];
  2033.     x2r = a[8] + a[24];
  2034.     x2i = a[9] + a[25];
  2035.     x3r = a[8] - a[24];
  2036.     x3i = a[9] - a[25];
  2037.     y0r = x0r + x2r;
  2038.     y0i = x0i + x2i;
  2039.     y4r = x0r - x2r;
  2040.     y4i = x0i - x2i;
  2041.     y8r = x1r - x3i;
  2042.     y8i = x1i + x3r;
  2043.     y12r = x1r + x3i;
  2044.     y12i = x1i - x3r;
  2045.     x0r = a[2] + a[18];
  2046.     x0i = a[3] + a[19];
  2047.     x1r = a[2] - a[18];
  2048.     x1i = a[3] - a[19];
  2049.     x2r = a[10] + a[26];
  2050.     x2i = a[11] + a[27];
  2051.     x3r = a[10] - a[26];
  2052.     x3i = a[11] - a[27];
  2053.     y1r = x0r + x2r;
  2054.     y1i = x0i + x2i;
  2055.     y5r = x0r - x2r;
  2056.     y5i = x0i - x2i;
  2057.     x0r = x1r - x3i;
  2058.     x0i = x1i + x3r;
  2059.     y9r = wk1r * x0r - wk1i * x0i;
  2060.     y9i = wk1r * x0i + wk1i * x0r;
  2061.     x0r = x1r + x3i;
  2062.     x0i = x1i - x3r;
  2063.     y13r = wk1i * x0r - wk1r * x0i;
  2064.     y13i = wk1i * x0i + wk1r * x0r;
  2065.     x0r = a[4] + a[20];
  2066.     x0i = a[5] + a[21];
  2067.     x1r = a[4] - a[20];
  2068.     x1i = a[5] - a[21];
  2069.     x2r = a[12] + a[28];
  2070.     x2i = a[13] + a[29];
  2071.     x3r = a[12] - a[28];
  2072.     x3i = a[13] - a[29];
  2073.     y2r = x0r + x2r;
  2074.     y2i = x0i + x2i;
  2075.     y6r = x0r - x2r;
  2076.     y6i = x0i - x2i;
  2077.     x0r = x1r - x3i;
  2078.     x0i = x1i + x3r;
  2079.     y10r = wn4r * (x0r - x0i);
  2080.     y10i = wn4r * (x0i + x0r);
  2081.     x0r = x1r + x3i;
  2082.     x0i = x1i - x3r;
  2083.     y14r = wn4r * (x0r + x0i);
  2084.     y14i = wn4r * (x0i - x0r);
  2085.     x0r = a[6] + a[22];
  2086.     x0i = a[7] + a[23];
  2087.     x1r = a[6] - a[22];
  2088.     x1i = a[7] - a[23];
  2089.     x2r = a[14] + a[30];
  2090.     x2i = a[15] + a[31];
  2091.     x3r = a[14] - a[30];
  2092.     x3i = a[15] - a[31];
  2093.     y3r = x0r + x2r;
  2094.     y3i = x0i + x2i;
  2095.     y7r = x0r - x2r;
  2096.     y7i = x0i - x2i;
  2097.     x0r = x1r - x3i;
  2098.     x0i = x1i + x3r;
  2099.     y11r = wk1i * x0r - wk1r * x0i;
  2100.     y11i = wk1i * x0i + wk1r * x0r;
  2101.     x0r = x1r + x3i;
  2102.     x0i = x1i - x3r;
  2103.     y15r = wk1r * x0r - wk1i * x0i;
  2104.     y15i = wk1r * x0i + wk1i * x0r;
  2105.     x0r = y12r - y14r;
  2106.     x0i = y12i - y14i;
  2107.     x1r = y12r + y14r;
  2108.     x1i = y12i + y14i;
  2109.     x2r = y13r - y15r;
  2110.     x2i = y13i - y15i;
  2111.     x3r = y13r + y15r;
  2112.     x3i = y13i + y15i;
  2113.     a[24] = x0r + x2r;
  2114.     a[25] = x0i + x2i;
  2115.     a[26] = x0r - x2r;
  2116.     a[27] = x0i - x2i;
  2117.     a[28] = x1r - x3i;
  2118.     a[29] = x1i + x3r;
  2119.     a[30] = x1r + x3i;
  2120.     a[31] = x1i - x3r;
  2121.     x0r = y8r + y10r;
  2122.     x0i = y8i + y10i;
  2123.     x1r = y8r - y10r;
  2124.     x1i = y8i - y10i;
  2125.     x2r = y9r + y11r;
  2126.     x2i = y9i + y11i;
  2127.     x3r = y9r - y11r;
  2128.     x3i = y9i - y11i;
  2129.     a[16] = x0r + x2r;
  2130.     a[17] = x0i + x2i;
  2131.     a[18] = x0r - x2r;
  2132.     a[19] = x0i - x2i;
  2133.     a[20] = x1r - x3i;
  2134.     a[21] = x1i + x3r;
  2135.     a[22] = x1r + x3i;
  2136.     a[23] = x1i - x3r;
  2137.     x0r = y5r - y7i;
  2138.     x0i = y5i + y7r;
  2139.     x2r = wn4r * (x0r - x0i);
  2140.     x2i = wn4r * (x0i + x0r);
  2141.     x0r = y5r + y7i;
  2142.     x0i = y5i - y7r;
  2143.     x3r = wn4r * (x0r - x0i);
  2144.     x3i = wn4r * (x0i + x0r);
  2145.     x0r = y4r - y6i;
  2146.     x0i = y4i + y6r;
  2147.     x1r = y4r + y6i;
  2148.     x1i = y4i - y6r;
  2149.     a[8] = x0r + x2r;
  2150.     a[9] = x0i + x2i;
  2151.     a[10] = x0r - x2r;
  2152.     a[11] = x0i - x2i;
  2153.     a[12] = x1r - x3i;
  2154.     a[13] = x1i + x3r;
  2155.     a[14] = x1r + x3i;
  2156.     a[15] = x1i - x3r;
  2157.     x0r = y0r + y2r;
  2158.     x0i = y0i + y2i;
  2159.     x1r = y0r - y2r;
  2160.     x1i = y0i - y2i;
  2161.     x2r = y1r + y3r;
  2162.     x2i = y1i + y3i;
  2163.     x3r = y1r - y3r;
  2164.     x3i = y1i - y3i;
  2165.     a[0] = x0r + x2r;
  2166.     a[1] = x0i + x2i;
  2167.     a[2] = x0r - x2r;
  2168.     a[3] = x0i - x2i;
  2169.     a[4] = x1r - x3i;
  2170.     a[5] = x1i + x3r;
  2171.     a[6] = x1r + x3i;
  2172.     a[7] = x1i - x3r;
  2173. }
  2174.  
  2175.  
  2176. void cftf162(double *a, double *w)
  2177. {
  2178.     double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
  2179.         x0r, x0i, x1r, x1i, x2r, x2i, 
  2180.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2181.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
  2182.         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
  2183.         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
  2184.     
  2185.     wn4r = w[1];
  2186.     wk1r = w[4];
  2187.     wk1i = w[5];
  2188.     wk3r = w[6];
  2189.     wk3i = w[7];
  2190.     wk2r = w[8];
  2191.     wk2i = w[9];
  2192.     x1r = a[0] - a[17];
  2193.     x1i = a[1] + a[16];
  2194.     x0r = a[8] - a[25];
  2195.     x0i = a[9] + a[24];
  2196.     x2r = wn4r * (x0r - x0i);
  2197.     x2i = wn4r * (x0i + x0r);
  2198.     y0r = x1r + x2r;
  2199.     y0i = x1i + x2i;
  2200.     y4r = x1r - x2r;
  2201.     y4i = x1i - x2i;
  2202.     x1r = a[0] + a[17];
  2203.     x1i = a[1] - a[16];
  2204.     x0r = a[8] + a[25];
  2205.     x0i = a[9] - a[24];
  2206.     x2r = wn4r * (x0r - x0i);
  2207.     x2i = wn4r * (x0i + x0r);
  2208.     y8r = x1r - x2i;
  2209.     y8i = x1i + x2r;
  2210.     y12r = x1r + x2i;
  2211.     y12i = x1i - x2r;
  2212.     x0r = a[2] - a[19];
  2213.     x0i = a[3] + a[18];
  2214.     x1r = wk1r * x0r - wk1i * x0i;
  2215.     x1i = wk1r * x0i + wk1i * x0r;
  2216.     x0r = a[10] - a[27];
  2217.     x0i = a[11] + a[26];
  2218.     x2r = wk3i * x0r - wk3r * x0i;
  2219.     x2i = wk3i * x0i + wk3r * x0r;
  2220.     y1r = x1r + x2r;
  2221.     y1i = x1i + x2i;
  2222.     y5r = x1r - x2r;
  2223.     y5i = x1i - x2i;
  2224.     x0r = a[2] + a[19];
  2225.     x0i = a[3] - a[18];
  2226.     x1r = wk3r * x0r - wk3i * x0i;
  2227.     x1i = wk3r * x0i + wk3i * x0r;
  2228.     x0r = a[10] + a[27];
  2229.     x0i = a[11] - a[26];
  2230.     x2r = wk1r * x0r + wk1i * x0i;
  2231.     x2i = wk1r * x0i - wk1i * x0r;
  2232.     y9r = x1r - x2r;
  2233.     y9i = x1i - x2i;
  2234.     y13r = x1r + x2r;
  2235.     y13i = x1i + x2i;
  2236.     x0r = a[4] - a[21];
  2237.     x0i = a[5] + a[20];
  2238.     x1r = wk2r * x0r - wk2i * x0i;
  2239.     x1i = wk2r * x0i + wk2i * x0r;
  2240.     x0r = a[12] - a[29];
  2241.     x0i = a[13] + a[28];
  2242.     x2r = wk2i * x0r - wk2r * x0i;
  2243.     x2i = wk2i * x0i + wk2r * x0r;
  2244.     y2r = x1r + x2r;
  2245.     y2i = x1i + x2i;
  2246.     y6r = x1r - x2r;
  2247.     y6i = x1i - x2i;
  2248.     x0r = a[4] + a[21];
  2249.     x0i = a[5] - a[20];
  2250.     x1r = wk2i * x0r - wk2r * x0i;
  2251.     x1i = wk2i * x0i + wk2r * x0r;
  2252.     x0r = a[12] + a[29];
  2253.     x0i = a[13] - a[28];
  2254.     x2r = wk2r * x0r - wk2i * x0i;
  2255.     x2i = wk2r * x0i + wk2i * x0r;
  2256.     y10r = x1r - x2r;
  2257.     y10i = x1i - x2i;
  2258.     y14r = x1r + x2r;
  2259.     y14i = x1i + x2i;
  2260.     x0r = a[6] - a[23];
  2261.     x0i = a[7] + a[22];
  2262.     x1r = wk3r * x0r - wk3i * x0i;
  2263.     x1i = wk3r * x0i + wk3i * x0r;
  2264.     x0r = a[14] - a[31];
  2265.     x0i = a[15] + a[30];
  2266.     x2r = wk1i * x0r - wk1r * x0i;
  2267.     x2i = wk1i * x0i + wk1r * x0r;
  2268.     y3r = x1r + x2r;
  2269.     y3i = x1i + x2i;
  2270.     y7r = x1r - x2r;
  2271.     y7i = x1i - x2i;
  2272.     x0r = a[6] + a[23];
  2273.     x0i = a[7] - a[22];
  2274.     x1r = wk1i * x0r + wk1r * x0i;
  2275.     x1i = wk1i * x0i - wk1r * x0r;
  2276.     x0r = a[14] + a[31];
  2277.     x0i = a[15] - a[30];
  2278.     x2r = wk3i * x0r - wk3r * x0i;
  2279.     x2i = wk3i * x0i + wk3r * x0r;
  2280.     y11r = x1r + x2r;
  2281.     y11i = x1i + x2i;
  2282.     y15r = x1r - x2r;
  2283.     y15i = x1i - x2i;
  2284.     x1r = y0r + y2r;
  2285.     x1i = y0i + y2i;
  2286.     x2r = y1r + y3r;
  2287.     x2i = y1i + y3i;
  2288.     a[0] = x1r + x2r;
  2289.     a[1] = x1i + x2i;
  2290.     a[2] = x1r - x2r;
  2291.     a[3] = x1i - x2i;
  2292.     x1r = y0r - y2r;
  2293.     x1i = y0i - y2i;
  2294.     x2r = y1r - y3r;
  2295.     x2i = y1i - y3i;
  2296.     a[4] = x1r - x2i;
  2297.     a[5] = x1i + x2r;
  2298.     a[6] = x1r + x2i;
  2299.     a[7] = x1i - x2r;
  2300.     x1r = y4r - y6i;
  2301.     x1i = y4i + y6r;
  2302.     x0r = y5r - y7i;
  2303.     x0i = y5i + y7r;
  2304.     x2r = wn4r * (x0r - x0i);
  2305.     x2i = wn4r * (x0i + x0r);
  2306.     a[8] = x1r + x2r;
  2307.     a[9] = x1i + x2i;
  2308.     a[10] = x1r - x2r;
  2309.     a[11] = x1i - x2i;
  2310.     x1r = y4r + y6i;
  2311.     x1i = y4i - y6r;
  2312.     x0r = y5r + y7i;
  2313.     x0i = y5i - y7r;
  2314.     x2r = wn4r * (x0r - x0i);
  2315.     x2i = wn4r * (x0i + x0r);
  2316.     a[12] = x1r - x2i;
  2317.     a[13] = x1i + x2r;
  2318.     a[14] = x1r + x2i;
  2319.     a[15] = x1i - x2r;
  2320.     x1r = y8r + y10r;
  2321.     x1i = y8i + y10i;
  2322.     x2r = y9r - y11r;
  2323.     x2i = y9i - y11i;
  2324.     a[16] = x1r + x2r;
  2325.     a[17] = x1i + x2i;
  2326.     a[18] = x1r - x2r;
  2327.     a[19] = x1i - x2i;
  2328.     x1r = y8r - y10r;
  2329.     x1i = y8i - y10i;
  2330.     x2r = y9r + y11r;
  2331.     x2i = y9i + y11i;
  2332.     a[20] = x1r - x2i;
  2333.     a[21] = x1i + x2r;
  2334.     a[22] = x1r + x2i;
  2335.     a[23] = x1i - x2r;
  2336.     x1r = y12r - y14i;
  2337.     x1i = y12i + y14r;
  2338.     x0r = y13r + y15i;
  2339.     x0i = y13i - y15r;
  2340.     x2r = wn4r * (x0r - x0i);
  2341.     x2i = wn4r * (x0i + x0r);
  2342.     a[24] = x1r + x2r;
  2343.     a[25] = x1i + x2i;
  2344.     a[26] = x1r - x2r;
  2345.     a[27] = x1i - x2i;
  2346.     x1r = y12r + y14i;
  2347.     x1i = y12i - y14r;
  2348.     x0r = y13r - y15i;
  2349.     x0i = y13i + y15r;
  2350.     x2r = wn4r * (x0r - x0i);
  2351.     x2i = wn4r * (x0i + x0r);
  2352.     a[28] = x1r - x2i;
  2353.     a[29] = x1i + x2r;
  2354.     a[30] = x1r + x2i;
  2355.     a[31] = x1i - x2r;
  2356. }
  2357.  
  2358.  
  2359. void cftf081(double *a, double *w)
  2360. {
  2361.     double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  2362.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2363.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
  2364.     
  2365.     wn4r = w[1];
  2366.     x0r = a[0] + a[8];
  2367.     x0i = a[1] + a[9];
  2368.     x1r = a[0] - a[8];
  2369.     x1i = a[1] - a[9];
  2370.     x2r = a[4] + a[12];
  2371.     x2i = a[5] + a[13];
  2372.     x3r = a[4] - a[12];
  2373.     x3i = a[5] - a[13];
  2374.     y0r = x0r + x2r;
  2375.     y0i = x0i + x2i;
  2376.     y2r = x0r - x2r;
  2377.     y2i = x0i - x2i;
  2378.     y1r = x1r - x3i;
  2379.     y1i = x1i + x3r;
  2380.     y3r = x1r + x3i;
  2381.     y3i = x1i - x3r;
  2382.     x0r = a[2] + a[10];
  2383.     x0i = a[3] + a[11];
  2384.     x1r = a[2] - a[10];
  2385.     x1i = a[3] - a[11];
  2386.     x2r = a[6] + a[14];
  2387.     x2i = a[7] + a[15];
  2388.     x3r = a[6] - a[14];
  2389.     x3i = a[7] - a[15];
  2390.     y4r = x0r + x2r;
  2391.     y4i = x0i + x2i;
  2392.     y6r = x0r - x2r;
  2393.     y6i = x0i - x2i;
  2394.     x0r = x1r - x3i;
  2395.     x0i = x1i + x3r;
  2396.     x2r = x1r + x3i;
  2397.     x2i = x1i - x3r;
  2398.     y5r = wn4r * (x0r - x0i);
  2399.     y5i = wn4r * (x0r + x0i);
  2400.     y7r = wn4r * (x2r - x2i);
  2401.     y7i = wn4r * (x2r + x2i);
  2402.     a[8] = y1r + y5r;
  2403.     a[9] = y1i + y5i;
  2404.     a[10] = y1r - y5r;
  2405.     a[11] = y1i - y5i;
  2406.     a[12] = y3r - y7i;
  2407.     a[13] = y3i + y7r;
  2408.     a[14] = y3r + y7i;
  2409.     a[15] = y3i - y7r;
  2410.     a[0] = y0r + y4r;
  2411.     a[1] = y0i + y4i;
  2412.     a[2] = y0r - y4r;
  2413.     a[3] = y0i - y4i;
  2414.     a[4] = y2r - y6i;
  2415.     a[5] = y2i + y6r;
  2416.     a[6] = y2r + y6i;
  2417.     a[7] = y2i - y6r;
  2418. }
  2419.  
  2420.  
  2421. void cftf082(double *a, double *w)
  2422. {
  2423.     double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
  2424.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2425.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
  2426.     
  2427.     wn4r = w[1];
  2428.     wk1r = w[4];
  2429.     wk1i = w[5];
  2430.     y0r = a[0] - a[9];
  2431.     y0i = a[1] + a[8];
  2432.     y1r = a[0] + a[9];
  2433.     y1i = a[1] - a[8];
  2434.     x0r = a[4] - a[13];
  2435.     x0i = a[5] + a[12];
  2436.     y2r = wn4r * (x0r - x0i);
  2437.     y2i = wn4r * (x0i + x0r);
  2438.     x0r = a[4] + a[13];
  2439.     x0i = a[5] - a[12];
  2440.     y3r = wn4r * (x0r - x0i);
  2441.     y3i = wn4r * (x0i + x0r);
  2442.     x0r = a[2] - a[11];
  2443.     x0i = a[3] + a[10];
  2444.     y4r = wk1r * x0r - wk1i * x0i;
  2445.     y4i = wk1r * x0i + wk1i * x0r;
  2446.     x0r = a[2] + a[11];
  2447.     x0i = a[3] - a[10];
  2448.     y5r = wk1i * x0r - wk1r * x0i;
  2449.     y5i = wk1i * x0i + wk1r * x0r;
  2450.     x0r = a[6] - a[15];
  2451.     x0i = a[7] + a[14];
  2452.     y6r = wk1i * x0r - wk1r * x0i;
  2453.     y6i = wk1i * x0i + wk1r * x0r;
  2454.     x0r = a[6] + a[15];
  2455.     x0i = a[7] - a[14];
  2456.     y7r = wk1r * x0r - wk1i * x0i;
  2457.     y7i = wk1r * x0i + wk1i * x0r;
  2458.     x0r = y0r + y2r;
  2459.     x0i = y0i + y2i;
  2460.     x1r = y4r + y6r;
  2461.     x1i = y4i + y6i;
  2462.     a[0] = x0r + x1r;
  2463.     a[1] = x0i + x1i;
  2464.     a[2] = x0r - x1r;
  2465.     a[3] = x0i - x1i;
  2466.     x0r = y0r - y2r;
  2467.     x0i = y0i - y2i;
  2468.     x1r = y4r - y6r;
  2469.     x1i = y4i - y6i;
  2470.     a[4] = x0r - x1i;
  2471.     a[5] = x0i + x1r;
  2472.     a[6] = x0r + x1i;
  2473.     a[7] = x0i - x1r;
  2474.     x0r = y1r - y3i;
  2475.     x0i = y1i + y3r;
  2476.     x1r = y5r - y7r;
  2477.     x1i = y5i - y7i;
  2478.     a[8] = x0r + x1r;
  2479.     a[9] = x0i + x1i;
  2480.     a[10] = x0r - x1r;
  2481.     a[11] = x0i - x1i;
  2482.     x0r = y1r + y3i;
  2483.     x0i = y1i - y3r;
  2484.     x1r = y5r + y7r;
  2485.     x1i = y5i + y7i;
  2486.     a[12] = x0r - x1i;
  2487.     a[13] = x0i + x1r;
  2488.     a[14] = x0r + x1i;
  2489.     a[15] = x0i - x1r;
  2490. }
  2491.  
  2492.  
  2493. void cftf040(double *a)
  2494. {
  2495.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  2496.     
  2497.     x0r = a[0] + a[4];
  2498.     x0i = a[1] + a[5];
  2499.     x1r = a[0] - a[4];
  2500.     x1i = a[1] - a[5];
  2501.     x2r = a[2] + a[6];
  2502.     x2i = a[3] + a[7];
  2503.     x3r = a[2] - a[6];
  2504.     x3i = a[3] - a[7];
  2505.     a[0] = x0r + x2r;
  2506.     a[1] = x0i + x2i;
  2507.     a[4] = x0r - x2r;
  2508.     a[5] = x0i - x2i;
  2509.     a[2] = x1r - x3i;
  2510.     a[3] = x1i + x3r;
  2511.     a[6] = x1r + x3i;
  2512.     a[7] = x1i - x3r;
  2513. }
  2514.  
  2515.  
  2516. void cftb040(double *a)
  2517. {
  2518.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  2519.     
  2520.     x0r = a[0] + a[4];
  2521.     x0i = a[1] + a[5];
  2522.     x1r = a[0] - a[4];
  2523.     x1i = a[1] - a[5];
  2524.     x2r = a[2] + a[6];
  2525.     x2i = a[3] + a[7];
  2526.     x3r = a[2] - a[6];
  2527.     x3i = a[3] - a[7];
  2528.     a[0] = x0r + x2r;
  2529.     a[1] = x0i + x2i;
  2530.     a[4] = x0r - x2r;
  2531.     a[5] = x0i - x2i;
  2532.     a[2] = x1r + x3i;
  2533.     a[3] = x1i - x3r;
  2534.     a[6] = x1r - x3i;
  2535.     a[7] = x1i + x3r;
  2536. }
  2537.  
  2538.  
  2539. void cftx020(double *a)
  2540. {
  2541.     double x0r, x0i;
  2542.     
  2543.     x0r = a[0] - a[2];
  2544.     x0i = a[1] - a[3];
  2545.     a[0] += a[2];
  2546.     a[1] += a[3];
  2547.     a[2] = x0r;
  2548.     a[3] = x0i;
  2549. }
  2550.  
  2551.  
  2552. void rftfsub(int n, double *a, int nc, double *c)
  2553. {
  2554.     int j, k, kk, ks, m;
  2555.     double wkr, wki, xr, xi, yr, yi;
  2556.     
  2557.     m = n >> 1;
  2558.     ks = 2 * nc / m;
  2559.     kk = 0;
  2560.     for (j = 2; j < m; j += 2) {
  2561.         k = n - j;
  2562.         kk += ks;
  2563.         wkr = 0.5 - c[nc - kk];
  2564.         wki = c[kk];
  2565.         xr = a[j] - a[k];
  2566.         xi = a[j + 1] + a[k + 1];
  2567.         yr = wkr * xr - wki * xi;
  2568.         yi = wkr * xi + wki * xr;
  2569.         a[j] -= yr;
  2570.         a[j + 1] -= yi;
  2571.         a[k] += yr;
  2572.         a[k + 1] -= yi;
  2573.     }
  2574. }
  2575.  
  2576.  
  2577. void rftbsub(int n, double *a, int nc, double *c)
  2578. {
  2579.     int j, k, kk, ks, m;
  2580.     double wkr, wki, xr, xi, yr, yi;
  2581.     
  2582.     m = n >> 1;
  2583.     ks = 2 * nc / m;
  2584.     kk = 0;
  2585.     for (j = 2; j < m; j += 2) {
  2586.         k = n - j;
  2587.         kk += ks;
  2588.         wkr = 0.5 - c[nc - kk];
  2589.         wki = c[kk];
  2590.         xr = a[j] - a[k];
  2591.         xi = a[j + 1] + a[k + 1];
  2592.         yr = wkr * xr + wki * xi;
  2593.         yi = wkr * xi - wki * xr;
  2594.         a[j] -= yr;
  2595.         a[j + 1] -= yi;
  2596.         a[k] += yr;
  2597.         a[k + 1] -= yi;
  2598.     }
  2599. }
  2600.  
  2601.  
  2602. void dctsub(int n, double *a, int nc, double *c)
  2603. {
  2604.     int j, k, kk, ks, m;
  2605.     double wkr, wki, xr;
  2606.     
  2607.     m = n >> 1;
  2608.     ks = nc / n;
  2609.     kk = 0;
  2610.     for (j = 1; j < m; j++) {
  2611.         k = n - j;
  2612.         kk += ks;
  2613.         wkr = c[kk] - c[nc - kk];
  2614.         wki = c[kk] + c[nc - kk];
  2615.         xr = wki * a[j] - wkr * a[k];
  2616.         a[j] = wkr * a[j] + wki * a[k];
  2617.         a[k] = xr;
  2618.     }
  2619.     a[m] *= c[0];
  2620. }
  2621.  
  2622.  
  2623. void dstsub(int n, double *a, int nc, double *c)
  2624. {
  2625.     int j, k, kk, ks, m;
  2626.     double wkr, wki, xr;
  2627.     
  2628.     m = n >> 1;
  2629.     ks = nc / n;
  2630.     kk = 0;
  2631.     for (j = 1; j < m; j++) {
  2632.         k = n - j;
  2633.         kk += ks;
  2634.         wkr = c[kk] - c[nc - kk];
  2635.         wki = c[kk] + c[nc - kk];
  2636.         xr = wki * a[k] - wkr * a[j];
  2637.         a[k] = wkr * a[k] + wki * a[j];
  2638.         a[j] = xr;
  2639.     }
  2640.     a[m] *= c[0];
  2641. }
  2642.  
  2643.