开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > midi > 查看源码
fftsg.c
资源名称：vlc-1.0.5.zip　[点击查看]
上传用户：kjfoods
上传日期：2020-07-06
资源大小：29949k
文件大小：88k
源码类别：

midi

开发平台：

Unix_Linux

fftsg.c：源码内容
  1. /*****************************************************************************
  2.  * fftsg.c:
  3.  *****************************************************************************
  4.  * Copyright (C) 2004 the VideoLAN team
  5.  * $Id: 7d6488f72d33db03df1b99522f68c951b239573c $
  6.  *
  7.  * Authors: Cyril Deguet <asmax@videolan.org>
  8.  *          code from projectM http://xmms-projectm.sourceforge.net
  9.  *
  10.  * This program is free software; you can redistribute it and/or modify
  11.  * it under the terms of the GNU General Public License as published by
  12.  * the Free Software Foundation; either version 2 of the License, or
  13.  * (at your option) any later version.
  14.  *
  15.  * This program is distributed in the hope that it will be useful,
  16.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18.  * GNU General Public License for more details.
  19.  *
  20.  * You should have received a copy of the GNU General Public License
  21.  * along with this program; if not, write to the Free Software
  22.  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston MA 02110-1301, USA.
  23.  *****************************************************************************/
  27. /*
  28. Fast Fourier/Cosine/Sine Transform
  29.     dimension   :one
  30.     data length :power of 2
  31.     decimation  :frequency
  32.     radix       :split-radix
  33.     data        :inplace
  34.     table       :use
  35. functions
  36.     cdft: Complex Discrete Fourier Transform
  37.     rdft: Real Discrete Fourier Transform
  38.     ddct: Discrete Cosine Transform
  39.     ddst: Discrete Sine Transform
  40.     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
  41.     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
  42. function prototypes
  43.     void cdft(int, int, double *, int *, double *);
  44.     void rdft(int, int, double *, int *, double *);
  45.     void ddct(int, int, double *, int *, double *);
  46.     void ddst(int, int, double *, int *, double *);
  47.     void dfct(int, double *, double *, int *, double *);
  48.     void dfst(int, double *, double *, int *, double *);
  49. macro definitions
  50.     USE_CDFT_PTHREADS : default=not defined
  51.         CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
  52.         CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
  53.     USE_CDFT_WINTHREADS : default=not defined
  54.         CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
  55.         CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288
  58. -------- Complex DFT (Discrete Fourier Transform) --------
  59.     [definition]
  60.         <case1>
  61.             X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
  62.         <case2>
  63.             X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
  64.         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
  65.     [usage]
  66.         <case1>
  67.             ip[0] = 0; // first time only
  68.             cdft(2*n, 1, a, ip, w);
  69.         <case2>
  70.             ip[0] = 0; // first time only
  71.             cdft(2*n, -1, a, ip, w);
  72.     [parameters]
  73.         2*n            :data length (int)
  74.                         n >= 1, n = power of 2
  75.         a[0...2*n-1]   :input/output data (double *)
  76.                         input data
  77.                             a[2*j] = Re(x[j]), 
  78.                             a[2*j+1] = Im(x[j]), 0<=j<n
  79.                         output data
  80.                             a[2*k] = Re(X[k]), 
  81.                             a[2*k+1] = Im(X[k]), 0<=k<n
  82.         ip[0...*]      :work area for bit reversal (int *)
  83.                         length of ip >= 2+sqrt(n)
  84.                         strictly, 
  85.                         length of ip >= 
  86.                             2+(1<<(int)(log(n+0.5)/log(2))/2).
  87.                         ip[0],ip[1] are pointers of the cos/sin table.
  88.         w[0...n/2-1]   :cos/sin table (double *)
  89.                         w[],ip[] are initialized if ip[0] == 0.
  90.     [remark]
  91.         Inverse of 
  92.             cdft(2*n, -1, a, ip, w);
  93.         is 
  94.             cdft(2*n, 1, a, ip, w);
  95.             for (j = 0; j <= 2 * n - 1; j++) {
  96.                 a[j] *= 1.0 / n;
  97.             }
  98.         .
  101. -------- Real DFT / Inverse of Real DFT --------
  102.     [definition]
  103.         <case1> RDFT
  104.             R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
  105.             I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
  106.         <case2> IRDFT (excluding scale)
  107.             a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
  108.                    sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
  109.                    sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
  110.     [usage]
  111.         <case1>
  112.             ip[0] = 0; // first time only
  113.             rdft(n, 1, a, ip, w);
  114.         <case2>
  115.             ip[0] = 0; // first time only
  116.             rdft(n, -1, a, ip, w);
  117.     [parameters]
  118.         n              :data length (int)
  119.                         n >= 2, n = power of 2
  120.         a[0...n-1]     :input/output data (double *)
  121.                         <case1>
  122.                             output data
  123.                                 a[2*k] = R[k], 0<=k<n/2
  124.                                 a[2*k+1] = I[k], 0<k<n/2
  125.                                 a[1] = R[n/2]
  126.                         <case2>
  127.                             input data
  128.                                 a[2*j] = R[j], 0<=j<n/2
  129.                                 a[2*j+1] = I[j], 0<j<n/2
  130.                                 a[1] = R[n/2]
  131.         ip[0...*]      :work area for bit reversal (int *)
  132.                         length of ip >= 2+sqrt(n/2)
  133.                         strictly, 
  134.                         length of ip >= 
  135.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  136.                         ip[0],ip[1] are pointers of the cos/sin table.
  137.         w[0...n/2-1]   :cos/sin table (double *)
  138.                         w[],ip[] are initialized if ip[0] == 0.
  139.     [remark]
  140.         Inverse of 
  141.             rdft(n, 1, a, ip, w);
  142.         is 
  143.             rdft(n, -1, a, ip, w);
  144.             for (j = 0; j <= n - 1; j++) {
  145.                 a[j] *= 2.0 / n;
  146.             }
  147.         .
  150. -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
  151.     [definition]
  152.         <case1> IDCT (excluding scale)
  153.             C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
  154.         <case2> DCT
  155.             C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
  156.     [usage]
  157.         <case1>
  158.             ip[0] = 0; // first time only
  159.             ddct(n, 1, a, ip, w);
  160.         <case2>
  161.             ip[0] = 0; // first time only
  162.             ddct(n, -1, a, ip, w);
  163.     [parameters]
  164.         n              :data length (int)
  165.                         n >= 2, n = power of 2
  166.         a[0...n-1]     :input/output data (double *)
  167.                         output data
  168.                             a[k] = C[k], 0<=k<n
  169.         ip[0...*]      :work area for bit reversal (int *)
  170.                         length of ip >= 2+sqrt(n/2)
  171.                         strictly, 
  172.                         length of ip >= 
  173.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  174.                         ip[0],ip[1] are pointers of the cos/sin table.
  175.         w[0...n*5/4-1] :cos/sin table (double *)
  176.                         w[],ip[] are initialized if ip[0] == 0.
  177.     [remark]
  178.         Inverse of 
  179.             ddct(n, -1, a, ip, w);
  180.         is 
  181.             a[0] *= 0.5;
  182.             ddct(n, 1, a, ip, w);
  183.             for (j = 0; j <= n - 1; j++) {
  184.                 a[j] *= 2.0 / n;
  185.             }
  186.         .
  189. -------- DST (Discrete Sine Transform) / Inverse of DST --------
  190.     [definition]
  191.         <case1> IDST (excluding scale)
  192.             S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
  193.         <case2> DST
  194.             S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
  195.     [usage]
  196.         <case1>
  197.             ip[0] = 0; // first time only
  198.             ddst(n, 1, a, ip, w);
  199.         <case2>
  200.             ip[0] = 0; // first time only
  201.             ddst(n, -1, a, ip, w);
  202.     [parameters]
  203.         n              :data length (int)
  204.                         n >= 2, n = power of 2
  205.         a[0...n-1]     :input/output data (double *)
  206.                         <case1>
  207.                             input data
  208.                                 a[j] = A[j], 0<j<n
  209.                                 a[0] = A[n]
  210.                             output data
  211.                                 a[k] = S[k], 0<=k<n
  212.                         <case2>
  213.                             output data
  214.                                 a[k] = S[k], 0<k<n
  215.                                 a[0] = S[n]
  216.         ip[0...*]      :work area for bit reversal (int *)
  217.                         length of ip >= 2+sqrt(n/2)
  218.                         strictly, 
  219.                         length of ip >= 
  220.                             2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  221.                         ip[0],ip[1] are pointers of the cos/sin table.
  222.         w[0...n*5/4-1] :cos/sin table (double *)
  223.                         w[],ip[] are initialized if ip[0] == 0.
  224.     [remark]
  225.         Inverse of 
  226.             ddst(n, -1, a, ip, w);
  227.         is 
  228.             a[0] *= 0.5;
  229.             ddst(n, 1, a, ip, w);
  230.             for (j = 0; j <= n - 1; j++) {
  231.                 a[j] *= 2.0 / n;
  232.             }
  233.         .
  236. -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
  237.     [definition]
  238.         C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
  239.     [usage]
  240.         ip[0] = 0; // first time only
  241.         dfct(n, a, t, ip, w);
  242.     [parameters]
  243.         n              :data length - 1 (int)
  244.                         n >= 2, n = power of 2
  245.         a[0...n]       :input/output data (double *)
  246.                         output data
  247.                             a[k] = C[k], 0<=k<=n
  248.         t[0...n/2]     :work area (double *)
  249.         ip[0...*]      :work area for bit reversal (int *)
  250.                         length of ip >= 2+sqrt(n/4)
  251.                         strictly, 
  252.                         length of ip >= 
  253.                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
  254.                         ip[0],ip[1] are pointers of the cos/sin table.
  255.         w[0...n*5/8-1] :cos/sin table (double *)
  256.                         w[],ip[] are initialized if ip[0] == 0.
  257.     [remark]
  258.         Inverse of 
  259.             a[0] *= 0.5;
  260.             a[n] *= 0.5;
  261.             dfct(n, a, t, ip, w);
  262.         is 
  263.             a[0] *= 0.5;
  264.             a[n] *= 0.5;
  265.             dfct(n, a, t, ip, w);
  266.             for (j = 0; j <= n; j++) {
  267.                 a[j] *= 2.0 / n;
  268.             }
  269.         .
  272. -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
  273.     [definition]
  274.         S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
  275.     [usage]
  276.         ip[0] = 0; // first time only
  277.         dfst(n, a, t, ip, w);
  278.     [parameters]
  279.         n              :data length + 1 (int)
  280.                         n >= 2, n = power of 2
  281.         a[0...n-1]     :input/output data (double *)
  282.                         output data
  283.                             a[k] = S[k], 0<k<n
  284.                         (a[0] is used for work area)
  285.         t[0...n/2-1]   :work area (double *)
  286.         ip[0...*]      :work area for bit reversal (int *)
  287.                         length of ip >= 2+sqrt(n/4)
  288.                         strictly, 
  289.                         length of ip >= 
  290.                             2+(1<<(int)(log(n/4+0.5)/log(2))/2).
  291.                         ip[0],ip[1] are pointers of the cos/sin table.
  292.         w[0...n*5/8-1] :cos/sin table (double *)
  293.                         w[],ip[] are initialized if ip[0] == 0.
  294.     [remark]
  295.         Inverse of 
  296.             dfst(n, a, t, ip, w);
  297.         is 
  298.             dfst(n, a, t, ip, w);
  299.             for (j = 1; j <= n - 1; j++) {
  300.                 a[j] *= 2.0 / n;
  301.             }
  302.         .
  305. Appendix :
  306.     The cos/sin table is recalculated when the larger table required.
  307.     w[] and ip[] are compatible with all routines.
  308. */
  311. void cdft(int n, int isgn, double *a, int *ip, double *w)
  312. {
  313.     void makewt(int nw, int *ip, double *w);
  314.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  315.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  316.     int nw;
  317.     
  318.     nw = ip[0];
  319.     if (n > (nw << 2)) {
  320.         nw = n >> 2;
  321.         makewt(nw, ip, w);
  322.     }
  323.     if (isgn >= 0) {
  324.         cftfsub(n, a, ip, nw, w);
  325.     } else {
  326.         cftbsub(n, a, ip, nw, w);
  327.     }
  328. }
  331. void rdft(int n, int isgn, double *a, int *ip, double *w)
  332. {
  333.     void makewt(int nw, int *ip, double *w);
  334.     void makect(int nc, int *ip, double *c);
  335.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  336.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  337.     void rftfsub(int n, double *a, int nc, double *c);
  338.     void rftbsub(int n, double *a, int nc, double *c);
  339.     int nw, nc;
  340.     double xi;
  341.     
  342.     nw = ip[0];
  343.     if (n > (nw << 2)) {
  344.         nw = n >> 2;
  345.         makewt(nw, ip, w);
  346.     }
  347.     nc = ip[1];
  348.     if (n > (nc << 2)) {
  349.         nc = n >> 2;
  350.         makect(nc, ip, w + nw);
  351.     }
  352.     if (isgn >= 0) {
  353.         if (n > 4) {
  354.             cftfsub(n, a, ip, nw, w);
  355.             rftfsub(n, a, nc, w + nw);
  356.         } else if (n == 4) {
  357.             cftfsub(n, a, ip, nw, w);
  358.         }
  359.         xi = a[0] - a[1];
  360.         a[0] += a[1];
  361.         a[1] = xi;
  362.     } else {
  363.         a[1] = 0.5 * (a[0] - a[1]);
  364.         a[0] -= a[1];
  365.         if (n > 4) {
  366.             rftbsub(n, a, nc, w + nw);
  367.             cftbsub(n, a, ip, nw, w);
  368.         } else if (n == 4) {
  369.             cftbsub(n, a, ip, nw, w);
  370.         }
  371.     }
  372. }
  375. void ddct(int n, int isgn, double *a, int *ip, double *w)
  376. {
  377.     void makewt(int nw, int *ip, double *w);
  378.     void makect(int nc, int *ip, double *c);
  379.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  380.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  381.     void rftfsub(int n, double *a, int nc, double *c);
  382.     void rftbsub(int n, double *a, int nc, double *c);
  383.     void dctsub(int n, double *a, int nc, double *c);
  384.     int j, nw, nc;
  385.     double xr;
  386.     
  387.     nw = ip[0];
  388.     if (n > (nw << 2)) {
  389.         nw = n >> 2;
  390.         makewt(nw, ip, w);
  391.     }
  392.     nc = ip[1];
  393.     if (n > nc) {
  394.         nc = n;
  395.         makect(nc, ip, w + nw);
  396.     }
  397.     if (isgn < 0) {
  398.         xr = a[n - 1];
  399.         for (j = n - 2; j >= 2; j -= 2) {
  400.             a[j + 1] = a[j] - a[j - 1];
  401.             a[j] += a[j - 1];
  402.         }
  403.         a[1] = a[0] - xr;
  404.         a[0] += xr;
  405.         if (n > 4) {
  406.             rftbsub(n, a, nc, w + nw);
  407.             cftbsub(n, a, ip, nw, w);
  408.         } else if (n == 4) {
  409.             cftbsub(n, a, ip, nw, w);
  410.         }
  411.     }
  412.     dctsub(n, a, nc, w + nw);
  413.     if (isgn >= 0) {
  414.         if (n > 4) {
  415.             cftfsub(n, a, ip, nw, w);
  416.             rftfsub(n, a, nc, w + nw);
  417.         } else if (n == 4) {
  418.             cftfsub(n, a, ip, nw, w);
  419.         }
  420.         xr = a[0] - a[1];
  421.         a[0] += a[1];
  422.         for (j = 2; j < n; j += 2) {
  423.             a[j - 1] = a[j] - a[j + 1];
  424.             a[j] += a[j + 1];
  425.         }
  426.         a[n - 1] = xr;
  427.     }
  428. }
  431. void ddst(int n, int isgn, double *a, int *ip, double *w)
  432. {
  433.     void makewt(int nw, int *ip, double *w);
  434.     void makect(int nc, int *ip, double *c);
  435.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  436.     void cftbsub(int n, double *a, int *ip, int nw, double *w);
  437.     void rftfsub(int n, double *a, int nc, double *c);
  438.     void rftbsub(int n, double *a, int nc, double *c);
  439.     void dstsub(int n, double *a, int nc, double *c);
  440.     int j, nw, nc;
  441.     double xr;
  442.     
  443.     nw = ip[0];
  444.     if (n > (nw << 2)) {
  445.         nw = n >> 2;
  446.         makewt(nw, ip, w);
  447.     }
  448.     nc = ip[1];
  449.     if (n > nc) {
  450.         nc = n;
  451.         makect(nc, ip, w + nw);
  452.     }
  453.     if (isgn < 0) {
  454.         xr = a[n - 1];
  455.         for (j = n - 2; j >= 2; j -= 2) {
  456.             a[j + 1] = -a[j] - a[j - 1];
  457.             a[j] -= a[j - 1];
  458.         }
  459.         a[1] = a[0] + xr;
  460.         a[0] -= xr;
  461.         if (n > 4) {
  462.             rftbsub(n, a, nc, w + nw);
  463.             cftbsub(n, a, ip, nw, w);
  464.         } else if (n == 4) {
  465.             cftbsub(n, a, ip, nw, w);
  466.         }
  467.     }
  468.     dstsub(n, a, nc, w + nw);
  469.     if (isgn >= 0) {
  470.         if (n > 4) {
  471.             cftfsub(n, a, ip, nw, w);
  472.             rftfsub(n, a, nc, w + nw);
  473.         } else if (n == 4) {
  474.             cftfsub(n, a, ip, nw, w);
  475.         }
  476.         xr = a[0] - a[1];
  477.         a[0] += a[1];
  478.         for (j = 2; j < n; j += 2) {
  479.             a[j - 1] = -a[j] - a[j + 1];
  480.             a[j] -= a[j + 1];
  481.         }
  482.         a[n - 1] = -xr;
  483.     }
  484. }
  487. void dfct(int n, double *a, double *t, int *ip, double *w)
  488. {
  489.     void makewt(int nw, int *ip, double *w);
  490.     void makect(int nc, int *ip, double *c);
  491.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  492.     void rftfsub(int n, double *a, int nc, double *c);
  493.     void dctsub(int n, double *a, int nc, double *c);
  494.     int j, k, l, m, mh, nw, nc;
  495.     double xr, xi, yr, yi;
  496.     
  497.     nw = ip[0];
  498.     if (n > (nw << 3)) {
  499.         nw = n >> 3;
  500.         makewt(nw, ip, w);
  501.     }
  502.     nc = ip[1];
  503.     if (n > (nc << 1)) {
  504.         nc = n >> 1;
  505.         makect(nc, ip, w + nw);
  506.     }
  507.     m = n >> 1;
  508.     yi = a[m];
  509.     xi = a[0] + a[n];
  510.     a[0] -= a[n];
  511.     t[0] = xi - yi;
  512.     t[m] = xi + yi;
  513.     if (n > 2) {
  514.         mh = m >> 1;
  515.         for (j = 1; j < mh; j++) {
  516.             k = m - j;
  517.             xr = a[j] - a[n - j];
  518.             xi = a[j] + a[n - j];
  519.             yr = a[k] - a[n - k];
  520.             yi = a[k] + a[n - k];
  521.             a[j] = xr;
  522.             a[k] = yr;
  523.             t[j] = xi - yi;
  524.             t[k] = xi + yi;
  525.         }
  526.         t[mh] = a[mh] + a[n - mh];
  527.         a[mh] -= a[n - mh];
  528.         dctsub(m, a, nc, w + nw);
  529.         if (m > 4) {
  530.             cftfsub(m, a, ip, nw, w);
  531.             rftfsub(m, a, nc, w + nw);
  532.         } else if (m == 4) {
  533.             cftfsub(m, a, ip, nw, w);
  534.         }
  535.         a[n - 1] = a[0] - a[1];
  536.         a[1] = a[0] + a[1];
  537.         for (j = m - 2; j >= 2; j -= 2) {
  538.             a[2 * j + 1] = a[j] + a[j + 1];
  539.             a[2 * j - 1] = a[j] - a[j + 1];
  540.         }
  541.         l = 2;
  542.         m = mh;
  543.         while (m >= 2) {
  544.             dctsub(m, t, nc, w + nw);
  545.             if (m > 4) {
  546.                 cftfsub(m, t, ip, nw, w);
  547.                 rftfsub(m, t, nc, w + nw);
  548.             } else if (m == 4) {
  549.                 cftfsub(m, t, ip, nw, w);
  550.             }
  551.             a[n - l] = t[0] - t[1];
  552.             a[l] = t[0] + t[1];
  553.             k = 0;
  554.             for (j = 2; j < m; j += 2) {
  555.                 k += l << 2;
  556.                 a[k - l] = t[j] - t[j + 1];
  557.                 a[k + l] = t[j] + t[j + 1];
  558.             }
  559.             l <<= 1;
  560.             mh = m >> 1;
  561.             for (j = 0; j < mh; j++) {
  562.                 k = m - j;
  563.                 t[j] = t[m + k] - t[m + j];
  564.                 t[k] = t[m + k] + t[m + j];
  565.             }
  566.             t[mh] = t[m + mh];
  567.             m = mh;
  568.         }
  569.         a[l] = t[0];
  570.         a[n] = t[2] - t[1];
  571.         a[0] = t[2] + t[1];
  572.     } else {
  573.         a[1] = a[0];
  574.         a[2] = t[0];
  575.         a[0] = t[1];
  576.     }
  577. }
  580. void dfst(int n, double *a, double *t, int *ip, double *w)
  581. {
  582.     void makewt(int nw, int *ip, double *w);
  583.     void makect(int nc, int *ip, double *c);
  584.     void cftfsub(int n, double *a, int *ip, int nw, double *w);
  585.     void rftfsub(int n, double *a, int nc, double *c);
  586.     void dstsub(int n, double *a, int nc, double *c);
  587.     int j, k, l, m, mh, nw, nc;
  588.     double xr, xi, yr, yi;
  589.     
  590.     nw = ip[0];
  591.     if (n > (nw << 3)) {
  592.         nw = n >> 3;
  593.         makewt(nw, ip, w);
  594.     }
  595.     nc = ip[1];
  596.     if (n > (nc << 1)) {
  597.         nc = n >> 1;
  598.         makect(nc, ip, w + nw);
  599.     }
  600.     if (n > 2) {
  601.         m = n >> 1;
  602.         mh = m >> 1;
  603.         for (j = 1; j < mh; j++) {
  604.             k = m - j;
  605.             xr = a[j] + a[n - j];
  606.             xi = a[j] - a[n - j];
  607.             yr = a[k] + a[n - k];
  608.             yi = a[k] - a[n - k];
  609.             a[j] = xr;
  610.             a[k] = yr;
  611.             t[j] = xi + yi;
  612.             t[k] = xi - yi;
  613.         }
  614.         t[0] = a[mh] - a[n - mh];
  615.         a[mh] += a[n - mh];
  616.         a[0] = a[m];
  617.         dstsub(m, a, nc, w + nw);
  618.         if (m > 4) {
  619.             cftfsub(m, a, ip, nw, w);
  620.             rftfsub(m, a, nc, w + nw);
  621.         } else if (m == 4) {
  622.             cftfsub(m, a, ip, nw, w);
  623.         }
  624.         a[n - 1] = a[1] - a[0];
  625.         a[1] = a[0] + a[1];
  626.         for (j = m - 2; j >= 2; j -= 2) {
  627.             a[2 * j + 1] = a[j] - a[j + 1];
  628.             a[2 * j - 1] = -a[j] - a[j + 1];
  629.         }
  630.         l = 2;
  631.         m = mh;
  632.         while (m >= 2) {
  633.             dstsub(m, t, nc, w + nw);
  634.             if (m > 4) {
  635.                 cftfsub(m, t, ip, nw, w);
  636.                 rftfsub(m, t, nc, w + nw);
  637.             } else if (m == 4) {
  638.                 cftfsub(m, t, ip, nw, w);
  639.             }
  640.             a[n - l] = t[1] - t[0];
  641.             a[l] = t[0] + t[1];
  642.             k = 0;
  643.             for (j = 2; j < m; j += 2) {
  644.                 k += l << 2;
  645.                 a[k - l] = -t[j] - t[j + 1];
  646.                 a[k + l] = t[j] - t[j + 1];
  647.             }
  648.             l <<= 1;
  649.             mh = m >> 1;
  650.             for (j = 1; j < mh; j++) {
  651.                 k = m - j;
  652.                 t[j] = t[m + k] + t[m + j];
  653.                 t[k] = t[m + k] - t[m + j];
  654.             }
  655.             t[0] = t[m + mh];
  656.             m = mh;
  657.         }
  658.         a[l] = t[0];
  659.     }
  660.     a[0] = 0;
  661. }
  664. /* -------- initializing routines -------- */
  667. #include <math.h>
  669. void makewt(int nw, int *ip, double *w)
  670. {
  671.     void makeipt(int nw, int *ip);
  672.     int j, nwh, nw0, nw1;
  673.     double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
  674.     
  675.     ip[0] = nw;
  676.     ip[1] = 1;
  677.     if (nw > 2) {
  678.         nwh = nw >> 1;
  679.         delta = atan(1.0) / nwh;
  680.         wn4r = cos(delta * nwh);
  681.         w[0] = 1;
  682.         w[1] = wn4r;
  683.         if (nwh == 4) {
  684.             w[2] = cos(delta * 2);
  685.             w[3] = sin(delta * 2);
  686.         } else if (nwh > 4) {
  687.             makeipt(nw, ip);
  688.             w[2] = 0.5 / cos(delta * 2);
  689.             w[3] = 0.5 / cos(delta * 6);
  690.             for (j = 4; j < nwh; j += 4) {
  691.                 w[j] = cos(delta * j);
  692.                 w[j + 1] = sin(delta * j);
  693.                 w[j + 2] = cos(3 * delta * j);
  694.                 w[j + 3] = -sin(3 * delta * j);
  695.             }
  696.         }
  697.         nw0 = 0;
  698.         while (nwh > 2) {
  699.             nw1 = nw0 + nwh;
  700.             nwh >>= 1;
  701.             w[nw1] = 1;
  702.             w[nw1 + 1] = wn4r;
  703.             if (nwh == 4) {
  704.                 wk1r = w[nw0 + 4];
  705.                 wk1i = w[nw0 + 5];
  706.                 w[nw1 + 2] = wk1r;
  707.                 w[nw1 + 3] = wk1i;
  708.             } else if (nwh > 4) {
  709.                 wk1r = w[nw0 + 4];
  710.                 wk3r = w[nw0 + 6];
  711.                 w[nw1 + 2] = 0.5 / wk1r;
  712.                 w[nw1 + 3] = 0.5 / wk3r;
  713.                 for (j = 4; j < nwh; j += 4) {
  714.                     wk1r = w[nw0 + 2 * j];
  715.                     wk1i = w[nw0 + 2 * j + 1];
  716.                     wk3r = w[nw0 + 2 * j + 2];
  717.                     wk3i = w[nw0 + 2 * j + 3];
  718.                     w[nw1 + j] = wk1r;
  719.                     w[nw1 + j + 1] = wk1i;
  720.                     w[nw1 + j + 2] = wk3r;
  721.                     w[nw1 + j + 3] = wk3i;
  722.                 }
  723.             }
  724.             nw0 = nw1;
  725.         }
  726.     }
  727. }
  730. void makeipt(int nw, int *ip)
  731. {
  732.     int j, l, m, m2, p, q;
  733.     
  734.     ip[2] = 0;
  735.     ip[3] = 16;
  736.     m = 2;
  737.     for (l = nw; l > 32; l >>= 2) {
  738.         m2 = m << 1;
  739.         q = m2 << 3;
  740.         for (j = m; j < m2; j++) {
  741.             p = ip[j] << 2;
  742.             ip[m + j] = p;
  743.             ip[m2 + j] = p + q;
  744.         }
  745.         m = m2;
  746.     }
  747. }
  750. void makect(int nc, int *ip, double *c)
  751. {
  752.     int j, nch;
  753.     double delta;
  754.     
  755.     ip[1] = nc;
  756.     if (nc > 1) {
  757.         nch = nc >> 1;
  758.         delta = atan(1.0) / nch;
  759.         c[0] = cos(delta * nch);
  760.         c[nch] = 0.5 * c[0];
  761.         for (j = 1; j < nch; j++) {
  762.             c[j] = 0.5 * cos(delta * j);
  763.             c[nc - j] = 0.5 * sin(delta * j);
  764.         }
  765.     }
  766. }
  769. /* -------- child routines -------- */
  772. #ifdef USE_CDFT_PTHREADS
  773. #define USE_CDFT_THREADS
  774. #ifndef CDFT_THREADS_BEGIN_N
  775. #define CDFT_THREADS_BEGIN_N 8192
  776. #endif
  777. #ifndef CDFT_4THREADS_BEGIN_N
  778. #define CDFT_4THREADS_BEGIN_N 65536
  779. #endif
  780. #include <pthread.h>
  781. #include <stdlib.h>
  782. #include <stdio.h>
  783. #define cdft_thread_t pthread_t
  784. #define cdft_thread_create(thp,func,argp) { 
  785.     if (pthread_create(thp, NULL, func, (void *) argp) != 0) { 
  786.         fprintf(stderr, "cdft thread errorn"); 
  787.         exit(1); 
  788.     } 
  789. }
  790. #define cdft_thread_wait(th) { 
  791.     if (pthread_join(th, NULL) != 0) { 
  792.         fprintf(stderr, "cdft thread errorn"); 
  793.         exit(1); 
  794.     } 
  795. }
  796. #endif /* USE_CDFT_PTHREADS */
  799. #ifdef USE_CDFT_WINTHREADS
  800. #define USE_CDFT_THREADS
  801. #ifndef CDFT_THREADS_BEGIN_N
  802. #define CDFT_THREADS_BEGIN_N 32768
  803. #endif
  804. #ifndef CDFT_4THREADS_BEGIN_N
  805. #define CDFT_4THREADS_BEGIN_N 524288
  806. #endif
  807. #include <windows.h>
  808. #include <stdlib.h>
  809. #include <stdio.h>
  810. #define cdft_thread_t HANDLE
  811. #define cdft_thread_create(thp,func,argp) { 
  812.     DWORD thid; 
  813.     *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); 
  814.     if (*(thp) == 0) { 
  815.         fprintf(stderr, "cdft thread errorn"); 
  816.         exit(1); 
  817.     } 
  818. }
  819. #define cdft_thread_wait(th) { 
  820.     WaitForSingleObject(th, INFINITE); 
  821.     CloseHandle(th); 
  822. }
  823. #endif /* USE_CDFT_WINTHREADS */
  826. void cftfsub(int n, double *a, int *ip, int nw, double *w)
  827. {
  828.     void bitrv2(int n, int *ip, double *a);
  829.     void bitrv216(double *a);
  830.     void bitrv208(double *a);
  831.     void cftf1st(int n, double *a, double *w);
  832.     void cftrec4(int n, double *a, int nw, double *w);
  833.     void cftleaf(int n, int isplt, double *a, int nw, double *w);
  834.     void cftfx41(int n, double *a, int nw, double *w);
  835.     void cftf161(double *a, double *w);
  836.     void cftf081(double *a, double *w);
  837.     void cftf040(double *a);
  838.     void cftx020(double *a);
  839. #ifdef USE_CDFT_THREADS
  840.     void cftrec4_th(int n, double *a, int nw, double *w);
  841. #endif /* USE_CDFT_THREADS */
  842.     
  843.     if (n > 8) {
  844.         if (n > 32) {
  845.             cftf1st(n, a, &w[nw - (n >> 2)]);
  846. #ifdef USE_CDFT_THREADS
  847.             if (n > CDFT_THREADS_BEGIN_N) {
  848.                 cftrec4_th(n, a, nw, w);
  849.             } else 
  850. #endif /* USE_CDFT_THREADS */
  851.             if (n > 512) {
  852.                 cftrec4(n, a, nw, w);
  853.             } else if (n > 128) {
  854.                 cftleaf(n, 1, a, nw, w);
  855.             } else {
  856.                 cftfx41(n, a, nw, w);
  857.             }
  858.             bitrv2(n, ip, a);
  859.         } else if (n == 32) {
  860.             cftf161(a, &w[nw - 8]);
  861.             bitrv216(a);
  862.         } else {
  863.             cftf081(a, w);
  864.             bitrv208(a);
  865.         }
  866.     } else if (n == 8) {
  867.         cftf040(a);
  868.     } else if (n == 4) {
  869.         cftx020(a);
  870.     }
  871. }
  874. void cftbsub(int n, double *a, int *ip, int nw, double *w)
  875. {
  876.     void bitrv2conj(int n, int *ip, double *a);
  877.     void bitrv216neg(double *a);
  878.     void bitrv208neg(double *a);
  879.     void cftb1st(int n, double *a, double *w);
  880.     void cftrec4(int n, double *a, int nw, double *w);
  881.     void cftleaf(int n, int isplt, double *a, int nw, double *w);
  882.     void cftfx41(int n, double *a, int nw, double *w);
  883.     void cftf161(double *a, double *w);
  884.     void cftf081(double *a, double *w);
  885.     void cftb040(double *a);
  886.     void cftx020(double *a);
  887. #ifdef USE_CDFT_THREADS
  888.     void cftrec4_th(int n, double *a, int nw, double *w);
  889. #endif /* USE_CDFT_THREADS */
  890.     
  891.     if (n > 8) {
  892.         if (n > 32) {
  893.             cftb1st(n, a, &w[nw - (n >> 2)]);
  894. #ifdef USE_CDFT_THREADS
  895.             if (n > CDFT_THREADS_BEGIN_N) {
  896.                 cftrec4_th(n, a, nw, w);
  897.             } else 
  898. #endif /* USE_CDFT_THREADS */
  899.             if (n > 512) {
  900.                 cftrec4(n, a, nw, w);
  901.             } else if (n > 128) {
  902.                 cftleaf(n, 1, a, nw, w);
  903.             } else {
  904.                 cftfx41(n, a, nw, w);
  905.             }
  906.             bitrv2conj(n, ip, a);
  907.         } else if (n == 32) {
  908.             cftf161(a, &w[nw - 8]);
  909.             bitrv216neg(a);
  910.         } else {
  911.             cftf081(a, w);
  912.             bitrv208neg(a);
  913.         }
  914.     } else if (n == 8) {
  915.         cftb040(a);
  916.     } else if (n == 4) {
  917.         cftx020(a);
  918.     }
  919. }
  922. void bitrv2(int n, int *ip, double *a)
  923. {
  924.     int j, j1, k, k1, l, m, nh, nm;
  925.     double xr, xi, yr, yi;
  926.     
  927.     m = 1;
  928.     for (l = n >> 2; l > 8; l >>= 2) {
  929.         m <<= 1;
  930.     }
  931.     nh = n >> 1;
  932.     nm = 4 * m;
  933.     if (l == 8) {
  934.         for (k = 0; k < m; k++) {
  935.             for (j = 0; j < k; j++) {
  936.                 j1 = 4 * j + 2 * ip[m + k];
  937.                 k1 = 4 * k + 2 * ip[m + j];
  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 += nm;
  947.                 k1 += 2 * nm;
  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 += nm;
  957.                 k1 -= nm;
  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.                 j1 += nm;
  967.                 k1 += 2 * nm;
  968.                 xr = a[j1];
  969.                 xi = a[j1 + 1];
  970.                 yr = a[k1];
  971.                 yi = a[k1 + 1];
  972.                 a[j1] = yr;
  973.                 a[j1 + 1] = yi;
  974.                 a[k1] = xr;
  975.                 a[k1 + 1] = xi;
  976.                 j1 += nh;
  977.                 k1 += 2;
  978.                 xr = a[j1];
  979.                 xi = a[j1 + 1];
  980.                 yr = a[k1];
  981.                 yi = a[k1 + 1];
  982.                 a[j1] = yr;
  983.                 a[j1 + 1] = yi;
  984.                 a[k1] = xr;
  985.                 a[k1 + 1] = xi;
  986.                 j1 -= nm;
  987.                 k1 -= 2 * nm;
  988.                 xr = a[j1];
  989.                 xi = a[j1 + 1];
  990.                 yr = a[k1];
  991.                 yi = a[k1 + 1];
  992.                 a[j1] = yr;
  993.                 a[j1 + 1] = yi;
  994.                 a[k1] = xr;
  995.                 a[k1 + 1] = xi;
  996.                 j1 -= nm;
  997.                 k1 += nm;
  998.                 xr = a[j1];
  999.                 xi = a[j1 + 1];
  1000.                 yr = a[k1];
  1001.                 yi = a[k1 + 1];
  1002.                 a[j1] = yr;
  1003.                 a[j1 + 1] = yi;
  1004.                 a[k1] = xr;
  1005.                 a[k1 + 1] = xi;
  1006.                 j1 -= nm;
  1007.                 k1 -= 2 * nm;
  1008.                 xr = a[j1];
  1009.                 xi = a[j1 + 1];
  1010.                 yr = a[k1];
  1011.                 yi = a[k1 + 1];
  1012.                 a[j1] = yr;
  1013.                 a[j1 + 1] = yi;
  1014.                 a[k1] = xr;
  1015.                 a[k1 + 1] = xi;
  1016.                 j1 += 2;
  1017.                 k1 += nh;
  1018.                 xr = a[j1];
  1019.                 xi = a[j1 + 1];
  1020.                 yr = a[k1];
  1021.                 yi = a[k1 + 1];
  1022.                 a[j1] = yr;
  1023.                 a[j1 + 1] = yi;
  1024.                 a[k1] = xr;
  1025.                 a[k1 + 1] = xi;
  1026.                 j1 += nm;
  1027.                 k1 += 2 * nm;
  1028.                 xr = a[j1];
  1029.                 xi = a[j1 + 1];
  1030.                 yr = a[k1];
  1031.                 yi = a[k1 + 1];
  1032.                 a[j1] = yr;
  1033.                 a[j1 + 1] = yi;
  1034.                 a[k1] = xr;
  1035.                 a[k1 + 1] = xi;
  1036.                 j1 += nm;
  1037.                 k1 -= nm;
  1038.                 xr = a[j1];
  1039.                 xi = a[j1 + 1];
  1040.                 yr = a[k1];
  1041.                 yi = a[k1 + 1];
  1042.                 a[j1] = yr;
  1043.                 a[j1 + 1] = yi;
  1044.                 a[k1] = xr;
  1045.                 a[k1 + 1] = xi;
  1046.                 j1 += nm;
  1047.                 k1 += 2 * nm;
  1048.                 xr = a[j1];
  1049.                 xi = a[j1 + 1];
  1050.                 yr = a[k1];
  1051.                 yi = a[k1 + 1];
  1052.                 a[j1] = yr;
  1053.                 a[j1 + 1] = yi;
  1054.                 a[k1] = xr;
  1055.                 a[k1 + 1] = xi;
  1056.                 j1 -= nh;
  1057.                 k1 -= 2;
  1058.                 xr = a[j1];
  1059.                 xi = a[j1 + 1];
  1060.                 yr = a[k1];
  1061.                 yi = a[k1 + 1];
  1062.                 a[j1] = yr;
  1063.                 a[j1 + 1] = yi;
  1064.                 a[k1] = xr;
  1065.                 a[k1 + 1] = xi;
  1066.                 j1 -= nm;
  1067.                 k1 -= 2 * nm;
  1068.                 xr = a[j1];
  1069.                 xi = a[j1 + 1];
  1070.                 yr = a[k1];
  1071.                 yi = a[k1 + 1];
  1072.                 a[j1] = yr;
  1073.                 a[j1 + 1] = yi;
  1074.                 a[k1] = xr;
  1075.                 a[k1 + 1] = xi;
  1076.                 j1 -= nm;
  1077.                 k1 += nm;
  1078.                 xr = a[j1];
  1079.                 xi = a[j1 + 1];
  1080.                 yr = a[k1];
  1081.                 yi = a[k1 + 1];
  1082.                 a[j1] = yr;
  1083.                 a[j1 + 1] = yi;
  1084.                 a[k1] = xr;
  1085.                 a[k1 + 1] = xi;
  1086.                 j1 -= nm;
  1087.                 k1 -= 2 * nm;
  1088.                 xr = a[j1];
  1089.                 xi = a[j1 + 1];
  1090.                 yr = a[k1];
  1091.                 yi = a[k1 + 1];
  1092.                 a[j1] = yr;
  1093.                 a[j1 + 1] = yi;
  1094.                 a[k1] = xr;
  1095.                 a[k1 + 1] = xi;
  1096.             }
  1097.             k1 = 4 * k + 2 * ip[m + k];
  1098.             j1 = k1 + 2;
  1099.             k1 += nh;
  1100.             xr = a[j1];
  1101.             xi = a[j1 + 1];
  1102.             yr = a[k1];
  1103.             yi = a[k1 + 1];
  1104.             a[j1] = yr;
  1105.             a[j1 + 1] = yi;
  1106.             a[k1] = xr;
  1107.             a[k1 + 1] = xi;
  1108.             j1 += nm;
  1109.             k1 += 2 * nm;
  1110.             xr = a[j1];
  1111.             xi = a[j1 + 1];
  1112.             yr = a[k1];
  1113.             yi = a[k1 + 1];
  1114.             a[j1] = yr;
  1115.             a[j1 + 1] = yi;
  1116.             a[k1] = xr;
  1117.             a[k1 + 1] = xi;
  1118.             j1 += nm;
  1119.             k1 -= nm;
  1120.             xr = a[j1];
  1121.             xi = a[j1 + 1];
  1122.             yr = a[k1];
  1123.             yi = a[k1 + 1];
  1124.             a[j1] = yr;
  1125.             a[j1 + 1] = yi;
  1126.             a[k1] = xr;
  1127.             a[k1 + 1] = xi;
  1128.             j1 -= 2;
  1129.             k1 -= nh;
  1130.             xr = a[j1];
  1131.             xi = a[j1 + 1];
  1132.             yr = a[k1];
  1133.             yi = a[k1 + 1];
  1134.             a[j1] = yr;
  1135.             a[j1 + 1] = yi;
  1136.             a[k1] = xr;
  1137.             a[k1 + 1] = xi;
  1138.             j1 += nh + 2;
  1139.             k1 += nh + 2;
  1140.             xr = a[j1];
  1141.             xi = a[j1 + 1];
  1142.             yr = a[k1];
  1143.             yi = a[k1 + 1];
  1144.             a[j1] = yr;
  1145.             a[j1 + 1] = yi;
  1146.             a[k1] = xr;
  1147.             a[k1 + 1] = xi;
  1148.             j1 -= nh - nm;
  1149.             k1 += 2 * nm - 2;
  1150.             xr = a[j1];
  1151.             xi = a[j1 + 1];
  1152.             yr = a[k1];
  1153.             yi = a[k1 + 1];
  1154.             a[j1] = yr;
  1155.             a[j1 + 1] = yi;
  1156.             a[k1] = xr;
  1157.             a[k1 + 1] = xi;
  1158.         }
  1159.     } else {
  1160.         for (k = 0; k < m; k++) {
  1161.             for (j = 0; j < k; j++) {
  1162.                 j1 = 4 * j + ip[m + k];
  1163.                 k1 = 4 * k + ip[m + j];
  1164.                 xr = a[j1];
  1165.                 xi = a[j1 + 1];
  1166.                 yr = a[k1];
  1167.                 yi = a[k1 + 1];
  1168.                 a[j1] = yr;
  1169.                 a[j1 + 1] = yi;
  1170.                 a[k1] = xr;
  1171.                 a[k1 + 1] = xi;
  1172.                 j1 += nm;
  1173.                 k1 += nm;
  1174.                 xr = a[j1];
  1175.                 xi = a[j1 + 1];
  1176.                 yr = a[k1];
  1177.                 yi = a[k1 + 1];
  1178.                 a[j1] = yr;
  1179.                 a[j1 + 1] = yi;
  1180.                 a[k1] = xr;
  1181.                 a[k1 + 1] = xi;
  1182.                 j1 += nh;
  1183.                 k1 += 2;
  1184.                 xr = a[j1];
  1185.                 xi = a[j1 + 1];
  1186.                 yr = a[k1];
  1187.                 yi = a[k1 + 1];
  1188.                 a[j1] = yr;
  1189.                 a[j1 + 1] = yi;
  1190.                 a[k1] = xr;
  1191.                 a[k1 + 1] = xi;
  1192.                 j1 -= nm;
  1193.                 k1 -= nm;
  1194.                 xr = a[j1];
  1195.                 xi = a[j1 + 1];
  1196.                 yr = a[k1];
  1197.                 yi = a[k1 + 1];
  1198.                 a[j1] = yr;
  1199.                 a[j1 + 1] = yi;
  1200.                 a[k1] = xr;
  1201.                 a[k1 + 1] = xi;
  1202.                 j1 += 2;
  1203.                 k1 += nh;
  1204.                 xr = a[j1];
  1205.                 xi = a[j1 + 1];
  1206.                 yr = a[k1];
  1207.                 yi = a[k1 + 1];
  1208.                 a[j1] = yr;
  1209.                 a[j1 + 1] = yi;
  1210.                 a[k1] = xr;
  1211.                 a[k1 + 1] = xi;
  1212.                 j1 += nm;
  1213.                 k1 += nm;
  1214.                 xr = a[j1];
  1215.                 xi = a[j1 + 1];
  1216.                 yr = a[k1];
  1217.                 yi = a[k1 + 1];
  1218.                 a[j1] = yr;
  1219.                 a[j1 + 1] = yi;
  1220.                 a[k1] = xr;
  1221.                 a[k1 + 1] = xi;
  1222.                 j1 -= nh;
  1223.                 k1 -= 2;
  1224.                 xr = a[j1];
  1225.                 xi = a[j1 + 1];
  1226.                 yr = a[k1];
  1227.                 yi = a[k1 + 1];
  1228.                 a[j1] = yr;
  1229.                 a[j1 + 1] = yi;
  1230.                 a[k1] = xr;
  1231.                 a[k1 + 1] = xi;
  1232.                 j1 -= nm;
  1233.                 k1 -= nm;
  1234.                 xr = a[j1];
  1235.                 xi = a[j1 + 1];
  1236.                 yr = a[k1];
  1237.                 yi = a[k1 + 1];
  1238.                 a[j1] = yr;
  1239.                 a[j1 + 1] = yi;
  1240.                 a[k1] = xr;
  1241.                 a[k1 + 1] = xi;
  1242.             }
  1243.             k1 = 4 * k + ip[m + k];
  1244.             j1 = k1 + 2;
  1245.             k1 += nh;
  1246.             xr = a[j1];
  1247.             xi = a[j1 + 1];
  1248.             yr = a[k1];
  1249.             yi = a[k1 + 1];
  1250.             a[j1] = yr;
  1251.             a[j1 + 1] = yi;
  1252.             a[k1] = xr;
  1253.             a[k1 + 1] = xi;
  1254.             j1 += nm;
  1255.             k1 += nm;
  1256.             xr = a[j1];
  1257.             xi = a[j1 + 1];
  1258.             yr = a[k1];
  1259.             yi = a[k1 + 1];
  1260.             a[j1] = yr;
  1261.             a[j1 + 1] = yi;
  1262.             a[k1] = xr;
  1263.             a[k1 + 1] = xi;
  1264.         }
  1265.     }
  1266. }
  1269. void bitrv2conj(int n, int *ip, double *a)
  1270. {
  1271.     int j, j1, k, k1, l, m, nh, nm;
  1272.     double xr, xi, yr, yi;
  1273.     
  1274.     m = 1;
  1275.     for (l = n >> 2; l > 8; l >>= 2) {
  1276.         m <<= 1;
  1277.     }
  1278.     nh = n >> 1;
  1279.     nm = 4 * m;
  1280.     if (l == 8) {
  1281.         for (k = 0; k < m; k++) {
  1282.             for (j = 0; j < k; j++) {
  1283.                 j1 = 4 * j + 2 * ip[m + k];
  1284.                 k1 = 4 * k + 2 * ip[m + j];
  1285.                 xr = a[j1];
  1286.                 xi = -a[j1 + 1];
  1287.                 yr = a[k1];
  1288.                 yi = -a[k1 + 1];
  1289.                 a[j1] = yr;
  1290.                 a[j1 + 1] = yi;
  1291.                 a[k1] = xr;
  1292.                 a[k1 + 1] = xi;
  1293.                 j1 += nm;
  1294.                 k1 += 2 * nm;
  1295.                 xr = a[j1];
  1296.                 xi = -a[j1 + 1];
  1297.                 yr = a[k1];
  1298.                 yi = -a[k1 + 1];
  1299.                 a[j1] = yr;
  1300.                 a[j1 + 1] = yi;
  1301.                 a[k1] = xr;
  1302.                 a[k1 + 1] = xi;
  1303.                 j1 += nm;
  1304.                 k1 -= nm;
  1305.                 xr = a[j1];
  1306.                 xi = -a[j1 + 1];
  1307.                 yr = a[k1];
  1308.                 yi = -a[k1 + 1];
  1309.                 a[j1] = yr;
  1310.                 a[j1 + 1] = yi;
  1311.                 a[k1] = xr;
  1312.                 a[k1 + 1] = xi;
  1313.                 j1 += nm;
  1314.                 k1 += 2 * nm;
  1315.                 xr = a[j1];
  1316.                 xi = -a[j1 + 1];
  1317.                 yr = a[k1];
  1318.                 yi = -a[k1 + 1];
  1319.                 a[j1] = yr;
  1320.                 a[j1 + 1] = yi;
  1321.                 a[k1] = xr;
  1322.                 a[k1 + 1] = xi;
  1323.                 j1 += nh;
  1324.                 k1 += 2;
  1325.                 xr = a[j1];
  1326.                 xi = -a[j1 + 1];
  1327.                 yr = a[k1];
  1328.                 yi = -a[k1 + 1];
  1329.                 a[j1] = yr;
  1330.                 a[j1 + 1] = yi;
  1331.                 a[k1] = xr;
  1332.                 a[k1 + 1] = xi;
  1333.                 j1 -= nm;
  1334.                 k1 -= 2 * nm;
  1335.                 xr = a[j1];
  1336.                 xi = -a[j1 + 1];
  1337.                 yr = a[k1];
  1338.                 yi = -a[k1 + 1];
  1339.                 a[j1] = yr;
  1340.                 a[j1 + 1] = yi;
  1341.                 a[k1] = xr;
  1342.                 a[k1 + 1] = xi;
  1343.                 j1 -= nm;
  1344.                 k1 += nm;
  1345.                 xr = a[j1];
  1346.                 xi = -a[j1 + 1];
  1347.                 yr = a[k1];
  1348.                 yi = -a[k1 + 1];
  1349.                 a[j1] = yr;
  1350.                 a[j1 + 1] = yi;
  1351.                 a[k1] = xr;
  1352.                 a[k1 + 1] = xi;
  1353.                 j1 -= nm;
  1354.                 k1 -= 2 * nm;
  1355.                 xr = a[j1];
  1356.                 xi = -a[j1 + 1];
  1357.                 yr = a[k1];
  1358.                 yi = -a[k1 + 1];
  1359.                 a[j1] = yr;
  1360.                 a[j1 + 1] = yi;
  1361.                 a[k1] = xr;
  1362.                 a[k1 + 1] = xi;
  1363.                 j1 += 2;
  1364.                 k1 += nh;
  1365.                 xr = a[j1];
  1366.                 xi = -a[j1 + 1];
  1367.                 yr = a[k1];
  1368.                 yi = -a[k1 + 1];
  1369.                 a[j1] = yr;
  1370.                 a[j1 + 1] = yi;
  1371.                 a[k1] = xr;
  1372.                 a[k1 + 1] = xi;
  1373.                 j1 += nm;
  1374.                 k1 += 2 * nm;
  1375.                 xr = a[j1];
  1376.                 xi = -a[j1 + 1];
  1377.                 yr = a[k1];
  1378.                 yi = -a[k1 + 1];
  1379.                 a[j1] = yr;
  1380.                 a[j1 + 1] = yi;
  1381.                 a[k1] = xr;
  1382.                 a[k1 + 1] = xi;
  1383.                 j1 += nm;
  1384.                 k1 -= nm;
  1385.                 xr = a[j1];
  1386.                 xi = -a[j1 + 1];
  1387.                 yr = a[k1];
  1388.                 yi = -a[k1 + 1];
  1389.                 a[j1] = yr;
  1390.                 a[j1 + 1] = yi;
  1391.                 a[k1] = xr;
  1392.                 a[k1 + 1] = xi;
  1393.                 j1 += nm;
  1394.                 k1 += 2 * nm;
  1395.                 xr = a[j1];
  1396.                 xi = -a[j1 + 1];
  1397.                 yr = a[k1];
  1398.                 yi = -a[k1 + 1];
  1399.                 a[j1] = yr;
  1400.                 a[j1 + 1] = yi;
  1401.                 a[k1] = xr;
  1402.                 a[k1 + 1] = xi;
  1403.                 j1 -= nh;
  1404.                 k1 -= 2;
  1405.                 xr = a[j1];
  1406.                 xi = -a[j1 + 1];
  1407.                 yr = a[k1];
  1408.                 yi = -a[k1 + 1];
  1409.                 a[j1] = yr;
  1410.                 a[j1 + 1] = yi;
  1411.                 a[k1] = xr;
  1412.                 a[k1 + 1] = xi;
  1413.                 j1 -= nm;
  1414.                 k1 -= 2 * nm;
  1415.                 xr = a[j1];
  1416.                 xi = -a[j1 + 1];
  1417.                 yr = a[k1];
  1418.                 yi = -a[k1 + 1];
  1419.                 a[j1] = yr;
  1420.                 a[j1 + 1] = yi;
  1421.                 a[k1] = xr;
  1422.                 a[k1 + 1] = xi;
  1423.                 j1 -= nm;
  1424.                 k1 += nm;
  1425.                 xr = a[j1];
  1426.                 xi = -a[j1 + 1];
  1427.                 yr = a[k1];
  1428.                 yi = -a[k1 + 1];
  1429.                 a[j1] = yr;
  1430.                 a[j1 + 1] = yi;
  1431.                 a[k1] = xr;
  1432.                 a[k1 + 1] = xi;
  1433.                 j1 -= nm;
  1434.                 k1 -= 2 * nm;
  1435.                 xr = a[j1];
  1436.                 xi = -a[j1 + 1];
  1437.                 yr = a[k1];
  1438.                 yi = -a[k1 + 1];
  1439.                 a[j1] = yr;
  1440.                 a[j1 + 1] = yi;
  1441.                 a[k1] = xr;
  1442.                 a[k1 + 1] = xi;
  1443.             }
  1444.             k1 = 4 * k + 2 * ip[m + k];
  1445.             j1 = k1 + 2;
  1446.             k1 += nh;
  1447.             a[j1 - 1] = -a[j1 - 1];
  1448.             xr = a[j1];
  1449.             xi = -a[j1 + 1];
  1450.             yr = a[k1];
  1451.             yi = -a[k1 + 1];
  1452.             a[j1] = yr;
  1453.             a[j1 + 1] = yi;
  1454.             a[k1] = xr;
  1455.             a[k1 + 1] = xi;
  1456.             a[k1 + 3] = -a[k1 + 3];
  1457.             j1 += nm;
  1458.             k1 += 2 * nm;
  1459.             xr = a[j1];
  1460.             xi = -a[j1 + 1];
  1461.             yr = a[k1];
  1462.             yi = -a[k1 + 1];
  1463.             a[j1] = yr;
  1464.             a[j1 + 1] = yi;
  1465.             a[k1] = xr;
  1466.             a[k1 + 1] = xi;
  1467.             j1 += nm;
  1468.             k1 -= nm;
  1469.             xr = a[j1];
  1470.             xi = -a[j1 + 1];
  1471.             yr = a[k1];
  1472.             yi = -a[k1 + 1];
  1473.             a[j1] = yr;
  1474.             a[j1 + 1] = yi;
  1475.             a[k1] = xr;
  1476.             a[k1 + 1] = xi;
  1477.             j1 -= 2;
  1478.             k1 -= nh;
  1479.             xr = a[j1];
  1480.             xi = -a[j1 + 1];
  1481.             yr = a[k1];
  1482.             yi = -a[k1 + 1];
  1483.             a[j1] = yr;
  1484.             a[j1 + 1] = yi;
  1485.             a[k1] = xr;
  1486.             a[k1 + 1] = xi;
  1487.             j1 += nh + 2;
  1488.             k1 += nh + 2;
  1489.             xr = a[j1];
  1490.             xi = -a[j1 + 1];
  1491.             yr = a[k1];
  1492.             yi = -a[k1 + 1];
  1493.             a[j1] = yr;
  1494.             a[j1 + 1] = yi;
  1495.             a[k1] = xr;
  1496.             a[k1 + 1] = xi;
  1497.             j1 -= nh - nm;
  1498.             k1 += 2 * nm - 2;
  1499.             a[j1 - 1] = -a[j1 - 1];
  1500.             xr = a[j1];
  1501.             xi = -a[j1 + 1];
  1502.             yr = a[k1];
  1503.             yi = -a[k1 + 1];
  1504.             a[j1] = yr;
  1505.             a[j1 + 1] = yi;
  1506.             a[k1] = xr;
  1507.             a[k1 + 1] = xi;
  1508.             a[k1 + 3] = -a[k1 + 3];
  1509.         }
  1510.     } else {
  1511.         for (k = 0; k < m; k++) {
  1512.             for (j = 0; j < k; j++) {
  1513.                 j1 = 4 * j + ip[m + k];
  1514.                 k1 = 4 * k + ip[m + j];
  1515.                 xr = a[j1];
  1516.                 xi = -a[j1 + 1];
  1517.                 yr = a[k1];
  1518.                 yi = -a[k1 + 1];
  1519.                 a[j1] = yr;
  1520.                 a[j1 + 1] = yi;
  1521.                 a[k1] = xr;
  1522.                 a[k1 + 1] = xi;
  1523.                 j1 += nm;
  1524.                 k1 += nm;
  1525.                 xr = a[j1];
  1526.                 xi = -a[j1 + 1];
  1527.                 yr = a[k1];
  1528.                 yi = -a[k1 + 1];
  1529.                 a[j1] = yr;
  1530.                 a[j1 + 1] = yi;
  1531.                 a[k1] = xr;
  1532.                 a[k1 + 1] = xi;
  1533.                 j1 += nh;
  1534.                 k1 += 2;
  1535.                 xr = a[j1];
  1536.                 xi = -a[j1 + 1];
  1537.                 yr = a[k1];
  1538.                 yi = -a[k1 + 1];
  1539.                 a[j1] = yr;
  1540.                 a[j1 + 1] = yi;
  1541.                 a[k1] = xr;
  1542.                 a[k1 + 1] = xi;
  1543.                 j1 -= nm;
  1544.                 k1 -= nm;
  1545.                 xr = a[j1];
  1546.                 xi = -a[j1 + 1];
  1547.                 yr = a[k1];
  1548.                 yi = -a[k1 + 1];
  1549.                 a[j1] = yr;
  1550.                 a[j1 + 1] = yi;
  1551.                 a[k1] = xr;
  1552.                 a[k1 + 1] = xi;
  1553.                 j1 += 2;
  1554.                 k1 += nh;
  1555.                 xr = a[j1];
  1556.                 xi = -a[j1 + 1];
  1557.                 yr = a[k1];
  1558.                 yi = -a[k1 + 1];
  1559.                 a[j1] = yr;
  1560.                 a[j1 + 1] = yi;
  1561.                 a[k1] = xr;
  1562.                 a[k1 + 1] = xi;
  1563.                 j1 += nm;
  1564.                 k1 += nm;
  1565.                 xr = a[j1];
  1566.                 xi = -a[j1 + 1];
  1567.                 yr = a[k1];
  1568.                 yi = -a[k1 + 1];
  1569.                 a[j1] = yr;
  1570.                 a[j1 + 1] = yi;
  1571.                 a[k1] = xr;
  1572.                 a[k1 + 1] = xi;
  1573.                 j1 -= nh;
  1574.                 k1 -= 2;
  1575.                 xr = a[j1];
  1576.                 xi = -a[j1 + 1];
  1577.                 yr = a[k1];
  1578.                 yi = -a[k1 + 1];
  1579.                 a[j1] = yr;
  1580.                 a[j1 + 1] = yi;
  1581.                 a[k1] = xr;
  1582.                 a[k1 + 1] = xi;
  1583.                 j1 -= nm;
  1584.                 k1 -= nm;
  1585.                 xr = a[j1];
  1586.                 xi = -a[j1 + 1];
  1587.                 yr = a[k1];
  1588.                 yi = -a[k1 + 1];
  1589.                 a[j1] = yr;
  1590.                 a[j1 + 1] = yi;
  1591.                 a[k1] = xr;
  1592.                 a[k1 + 1] = xi;
  1593.             }
  1594.             k1 = 4 * k + ip[m + k];
  1595.             j1 = k1 + 2;
  1596.             k1 += nh;
  1597.             a[j1 - 1] = -a[j1 - 1];
  1598.             xr = a[j1];
  1599.             xi = -a[j1 + 1];
  1600.             yr = a[k1];
  1601.             yi = -a[k1 + 1];
  1602.             a[j1] = yr;
  1603.             a[j1 + 1] = yi;
  1604.             a[k1] = xr;
  1605.             a[k1 + 1] = xi;
  1606.             a[k1 + 3] = -a[k1 + 3];
  1607.             j1 += nm;
  1608.             k1 += nm;
  1609.             a[j1 - 1] = -a[j1 - 1];
  1610.             xr = a[j1];
  1611.             xi = -a[j1 + 1];
  1612.             yr = a[k1];
  1613.             yi = -a[k1 + 1];
  1614.             a[j1] = yr;
  1615.             a[j1 + 1] = yi;
  1616.             a[k1] = xr;
  1617.             a[k1 + 1] = xi;
  1618.             a[k1 + 3] = -a[k1 + 3];
  1619.         }
  1620.     }
  1621. }
  1624. void bitrv216(double *a)
  1625. {
  1626.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1627.         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
  1628.         x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
  1629.     
  1630.     x1r = a[2];
  1631.     x1i = a[3];
  1632.     x2r = a[4];
  1633.     x2i = a[5];
  1634.     x3r = a[6];
  1635.     x3i = a[7];
  1636.     x4r = a[8];
  1637.     x4i = a[9];
  1638.     x5r = a[10];
  1639.     x5i = a[11];
  1640.     x7r = a[14];
  1641.     x7i = a[15];
  1642.     x8r = a[16];
  1643.     x8i = a[17];
  1644.     x10r = a[20];
  1645.     x10i = a[21];
  1646.     x11r = a[22];
  1647.     x11i = a[23];
  1648.     x12r = a[24];
  1649.     x12i = a[25];
  1650.     x13r = a[26];
  1651.     x13i = a[27];
  1652.     x14r = a[28];
  1653.     x14i = a[29];
  1654.     a[2] = x8r;
  1655.     a[3] = x8i;
  1656.     a[4] = x4r;
  1657.     a[5] = x4i;
  1658.     a[6] = x12r;
  1659.     a[7] = x12i;
  1660.     a[8] = x2r;
  1661.     a[9] = x2i;
  1662.     a[10] = x10r;
  1663.     a[11] = x10i;
  1664.     a[14] = x14r;
  1665.     a[15] = x14i;
  1666.     a[16] = x1r;
  1667.     a[17] = x1i;
  1668.     a[20] = x5r;
  1669.     a[21] = x5i;
  1670.     a[22] = x13r;
  1671.     a[23] = x13i;
  1672.     a[24] = x3r;
  1673.     a[25] = x3i;
  1674.     a[26] = x11r;
  1675.     a[27] = x11i;
  1676.     a[28] = x7r;
  1677.     a[29] = x7i;
  1678. }
  1681. void bitrv216neg(double *a)
  1682. {
  1683.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1684.         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
  1685.         x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i, 
  1686.         x13r, x13i, x14r, x14i, x15r, x15i;
  1687.     
  1688.     x1r = a[2];
  1689.     x1i = a[3];
  1690.     x2r = a[4];
  1691.     x2i = a[5];
  1692.     x3r = a[6];
  1693.     x3i = a[7];
  1694.     x4r = a[8];
  1695.     x4i = a[9];
  1696.     x5r = a[10];
  1697.     x5i = a[11];
  1698.     x6r = a[12];
  1699.     x6i = a[13];
  1700.     x7r = a[14];
  1701.     x7i = a[15];
  1702.     x8r = a[16];
  1703.     x8i = a[17];
  1704.     x9r = a[18];
  1705.     x9i = a[19];
  1706.     x10r = a[20];
  1707.     x10i = a[21];
  1708.     x11r = a[22];
  1709.     x11i = a[23];
  1710.     x12r = a[24];
  1711.     x12i = a[25];
  1712.     x13r = a[26];
  1713.     x13i = a[27];
  1714.     x14r = a[28];
  1715.     x14i = a[29];
  1716.     x15r = a[30];
  1717.     x15i = a[31];
  1718.     a[2] = x15r;
  1719.     a[3] = x15i;
  1720.     a[4] = x7r;
  1721.     a[5] = x7i;
  1722.     a[6] = x11r;
  1723.     a[7] = x11i;
  1724.     a[8] = x3r;
  1725.     a[9] = x3i;
  1726.     a[10] = x13r;
  1727.     a[11] = x13i;
  1728.     a[12] = x5r;
  1729.     a[13] = x5i;
  1730.     a[14] = x9r;
  1731.     a[15] = x9i;
  1732.     a[16] = x1r;
  1733.     a[17] = x1i;
  1734.     a[18] = x14r;
  1735.     a[19] = x14i;
  1736.     a[20] = x6r;
  1737.     a[21] = x6i;
  1738.     a[22] = x10r;
  1739.     a[23] = x10i;
  1740.     a[24] = x2r;
  1741.     a[25] = x2i;
  1742.     a[26] = x12r;
  1743.     a[27] = x12i;
  1744.     a[28] = x4r;
  1745.     a[29] = x4i;
  1746.     a[30] = x8r;
  1747.     a[31] = x8i;
  1748. }
  1751. void bitrv208(double *a)
  1752. {
  1753.     double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
  1754.     
  1755.     x1r = a[2];
  1756.     x1i = a[3];
  1757.     x3r = a[6];
  1758.     x3i = a[7];
  1759.     x4r = a[8];
  1760.     x4i = a[9];
  1761.     x6r = a[12];
  1762.     x6i = a[13];
  1763.     a[2] = x4r;
  1764.     a[3] = x4i;
  1765.     a[6] = x6r;
  1766.     a[7] = x6i;
  1767.     a[8] = x1r;
  1768.     a[9] = x1i;
  1769.     a[12] = x3r;
  1770.     a[13] = x3i;
  1771. }
  1774. void bitrv208neg(double *a)
  1775. {
  1776.     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
  1777.         x5r, x5i, x6r, x6i, x7r, x7i;
  1778.     
  1779.     x1r = a[2];
  1780.     x1i = a[3];
  1781.     x2r = a[4];
  1782.     x2i = a[5];
  1783.     x3r = a[6];
  1784.     x3i = a[7];
  1785.     x4r = a[8];
  1786.     x4i = a[9];
  1787.     x5r = a[10];
  1788.     x5i = a[11];
  1789.     x6r = a[12];
  1790.     x6i = a[13];
  1791.     x7r = a[14];
  1792.     x7i = a[15];
  1793.     a[2] = x7r;
  1794.     a[3] = x7i;
  1795.     a[4] = x3r;
  1796.     a[5] = x3i;
  1797.     a[6] = x5r;
  1798.     a[7] = x5i;
  1799.     a[8] = x1r;
  1800.     a[9] = x1i;
  1801.     a[10] = x6r;
  1802.     a[11] = x6i;
  1803.     a[12] = x2r;
  1804.     a[13] = x2i;
  1805.     a[14] = x4r;
  1806.     a[15] = x4i;
  1807. }
  1810. void cftf1st(int n, double *a, double *w)
  1811. {
  1812.     int j, j0, j1, j2, j3, k, m, mh;
  1813.     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
  1814.         wd1r, wd1i, wd3r, wd3i;
  1815.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  1816.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
  1817.     
  1818.     mh = n >> 3;
  1819.     m = 2 * mh;
  1820.     j1 = m;
  1821.     j2 = j1 + m;
  1822.     j3 = j2 + m;
  1823.     x0r = a[0] + a[j2];
  1824.     x0i = a[1] + a[j2 + 1];
  1825.     x1r = a[0] - a[j2];
  1826.     x1i = a[1] - a[j2 + 1];
  1827.     x2r = a[j1] + a[j3];
  1828.     x2i = a[j1 + 1] + a[j3 + 1];
  1829.     x3r = a[j1] - a[j3];
  1830.     x3i = a[j1 + 1] - a[j3 + 1];
  1831.     a[0] = x0r + x2r;
  1832.     a[1] = x0i + x2i;
  1833.     a[j1] = x0r - x2r;
  1834.     a[j1 + 1] = x0i - x2i;
  1835.     a[j2] = x1r - x3i;
  1836.     a[j2 + 1] = x1i + x3r;
  1837.     a[j3] = x1r + x3i;
  1838.     a[j3 + 1] = x1i - x3r;
  1839.     wn4r = w[1];
  1840.     csc1 = w[2];
  1841.     csc3 = w[3];
  1842.     wd1r = 1;
  1843.     wd1i = 0;
  1844.     wd3r = 1;
  1845.     wd3i = 0;
  1846.     k = 0;
  1847.     for (j = 2; j < mh - 2; j += 4) {
  1848.         k += 4;
  1849.         wk1r = csc1 * (wd1r + w[k]);
  1850.         wk1i = csc1 * (wd1i + w[k + 1]);
  1851.         wk3r = csc3 * (wd3r + w[k + 2]);
  1852.         wk3i = csc3 * (wd3i + w[k + 3]);
  1853.         wd1r = w[k];
  1854.         wd1i = w[k + 1];
  1855.         wd3r = w[k + 2];
  1856.         wd3i = w[k + 3];
  1857.         j1 = j + m;
  1858.         j2 = j1 + m;
  1859.         j3 = j2 + m;
  1860.         x0r = a[j] + a[j2];
  1861.         x0i = a[j + 1] + a[j2 + 1];
  1862.         x1r = a[j] - a[j2];
  1863.         x1i = a[j + 1] - a[j2 + 1];
  1864.         y0r = a[j + 2] + a[j2 + 2];
  1865.         y0i = a[j + 3] + a[j2 + 3];
  1866.         y1r = a[j + 2] - a[j2 + 2];
  1867.         y1i = a[j + 3] - a[j2 + 3];
  1868.         x2r = a[j1] + a[j3];
  1869.         x2i = a[j1 + 1] + a[j3 + 1];
  1870.         x3r = a[j1] - a[j3];
  1871.         x3i = a[j1 + 1] - a[j3 + 1];
  1872.         y2r = a[j1 + 2] + a[j3 + 2];
  1873.         y2i = a[j1 + 3] + a[j3 + 3];
  1874.         y3r = a[j1 + 2] - a[j3 + 2];
  1875.         y3i = a[j1 + 3] - a[j3 + 3];
  1876.         a[j] = x0r + x2r;
  1877.         a[j + 1] = x0i + x2i;
  1878.         a[j + 2] = y0r + y2r;
  1879.         a[j + 3] = y0i + y2i;
  1880.         a[j1] = x0r - x2r;
  1881.         a[j1 + 1] = x0i - x2i;
  1882.         a[j1 + 2] = y0r - y2r;
  1883.         a[j1 + 3] = y0i - y2i;
  1884.         x0r = x1r - x3i;
  1885.         x0i = x1i + x3r;
  1886.         a[j2] = wk1r * x0r - wk1i * x0i;
  1887.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  1888.         x0r = y1r - y3i;
  1889.         x0i = y1i + y3r;
  1890.         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
  1891.         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
  1892.         x0r = x1r + x3i;
  1893.         x0i = x1i - x3r;
  1894.         a[j3] = wk3r * x0r + wk3i * x0i;
  1895.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  1896.         x0r = y1r + y3i;
  1897.         x0i = y1i - y3r;
  1898.         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
  1899.         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
  1900.         j0 = m - j;
  1901.         j1 = j0 + m;
  1902.         j2 = j1 + m;
  1903.         j3 = j2 + m;
  1904.         x0r = a[j0] + a[j2];
  1905.         x0i = a[j0 + 1] + a[j2 + 1];
  1906.         x1r = a[j0] - a[j2];
  1907.         x1i = a[j0 + 1] - a[j2 + 1];
  1908.         y0r = a[j0 - 2] + a[j2 - 2];
  1909.         y0i = a[j0 - 1] + a[j2 - 1];
  1910.         y1r = a[j0 - 2] - a[j2 - 2];
  1911.         y1i = a[j0 - 1] - a[j2 - 1];
  1912.         x2r = a[j1] + a[j3];
  1913.         x2i = a[j1 + 1] + a[j3 + 1];
  1914.         x3r = a[j1] - a[j3];
  1915.         x3i = a[j1 + 1] - a[j3 + 1];
  1916.         y2r = a[j1 - 2] + a[j3 - 2];
  1917.         y2i = a[j1 - 1] + a[j3 - 1];
  1918.         y3r = a[j1 - 2] - a[j3 - 2];
  1919.         y3i = a[j1 - 1] - a[j3 - 1];
  1920.         a[j0] = x0r + x2r;
  1921.         a[j0 + 1] = x0i + x2i;
  1922.         a[j0 - 2] = y0r + y2r;
  1923.         a[j0 - 1] = y0i + y2i;
  1924.         a[j1] = x0r - x2r;
  1925.         a[j1 + 1] = x0i - x2i;
  1926.         a[j1 - 2] = y0r - y2r;
  1927.         a[j1 - 1] = y0i - y2i;
  1928.         x0r = x1r - x3i;
  1929.         x0i = x1i + x3r;
  1930.         a[j2] = wk1i * x0r - wk1r * x0i;
  1931.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  1932.         x0r = y1r - y3i;
  1933.         x0i = y1i + y3r;
  1934.         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
  1935.         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
  1936.         x0r = x1r + x3i;
  1937.         x0i = x1i - x3r;
  1938.         a[j3] = wk3i * x0r + wk3r * x0i;
  1939.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  1940.         x0r = y1r + y3i;
  1941.         x0i = y1i - y3r;
  1942.         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
  1943.         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
  1944.     }
  1945.     wk1r = csc1 * (wd1r + wn4r);
  1946.     wk1i = csc1 * (wd1i + wn4r);
  1947.     wk3r = csc3 * (wd3r - wn4r);
  1948.     wk3i = csc3 * (wd3i - wn4r);
  1949.     j0 = mh;
  1950.     j1 = j0 + m;
  1951.     j2 = j1 + m;
  1952.     j3 = j2 + m;
  1953.     x0r = a[j0 - 2] + a[j2 - 2];
  1954.     x0i = a[j0 - 1] + a[j2 - 1];
  1955.     x1r = a[j0 - 2] - a[j2 - 2];
  1956.     x1i = a[j0 - 1] - a[j2 - 1];
  1957.     x2r = a[j1 - 2] + a[j3 - 2];
  1958.     x2i = a[j1 - 1] + a[j3 - 1];
  1959.     x3r = a[j1 - 2] - a[j3 - 2];
  1960.     x3i = a[j1 - 1] - a[j3 - 1];
  1961.     a[j0 - 2] = x0r + x2r;
  1962.     a[j0 - 1] = x0i + x2i;
  1963.     a[j1 - 2] = x0r - x2r;
  1964.     a[j1 - 1] = x0i - x2i;
  1965.     x0r = x1r - x3i;
  1966.     x0i = x1i + x3r;
  1967.     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
  1968.     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
  1969.     x0r = x1r + x3i;
  1970.     x0i = x1i - x3r;
  1971.     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
  1972.     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
  1973.     x0r = a[j0] + a[j2];
  1974.     x0i = a[j0 + 1] + a[j2 + 1];
  1975.     x1r = a[j0] - a[j2];
  1976.     x1i = a[j0 + 1] - a[j2 + 1];
  1977.     x2r = a[j1] + a[j3];
  1978.     x2i = a[j1 + 1] + a[j3 + 1];
  1979.     x3r = a[j1] - a[j3];
  1980.     x3i = a[j1 + 1] - a[j3 + 1];
  1981.     a[j0] = x0r + x2r;
  1982.     a[j0 + 1] = x0i + x2i;
  1983.     a[j1] = x0r - x2r;
  1984.     a[j1 + 1] = x0i - x2i;
  1985.     x0r = x1r - x3i;
  1986.     x0i = x1i + x3r;
  1987.     a[j2] = wn4r * (x0r - x0i);
  1988.     a[j2 + 1] = wn4r * (x0i + x0r);
  1989.     x0r = x1r + x3i;
  1990.     x0i = x1i - x3r;
  1991.     a[j3] = -wn4r * (x0r + x0i);
  1992.     a[j3 + 1] = -wn4r * (x0i - x0r);
  1993.     x0r = a[j0 + 2] + a[j2 + 2];
  1994.     x0i = a[j0 + 3] + a[j2 + 3];
  1995.     x1r = a[j0 + 2] - a[j2 + 2];
  1996.     x1i = a[j0 + 3] - a[j2 + 3];
  1997.     x2r = a[j1 + 2] + a[j3 + 2];
  1998.     x2i = a[j1 + 3] + a[j3 + 3];
  1999.     x3r = a[j1 + 2] - a[j3 + 2];
  2000.     x3i = a[j1 + 3] - a[j3 + 3];
  2001.     a[j0 + 2] = x0r + x2r;
  2002.     a[j0 + 3] = x0i + x2i;
  2003.     a[j1 + 2] = x0r - x2r;
  2004.     a[j1 + 3] = x0i - x2i;
  2005.     x0r = x1r - x3i;
  2006.     x0i = x1i + x3r;
  2007.     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
  2008.     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
  2009.     x0r = x1r + x3i;
  2010.     x0i = x1i - x3r;
  2011.     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
  2012.     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
  2013. }
  2016. void cftb1st(int n, double *a, double *w)
  2017. {
  2018.     int j, j0, j1, j2, j3, k, m, mh;
  2019.     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
  2020.         wd1r, wd1i, wd3r, wd3i;
  2021.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  2022.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
  2023.     
  2024.     mh = n >> 3;
  2025.     m = 2 * mh;
  2026.     j1 = m;
  2027.     j2 = j1 + m;
  2028.     j3 = j2 + m;
  2029.     x0r = a[0] + a[j2];
  2030.     x0i = -a[1] - a[j2 + 1];
  2031.     x1r = a[0] - a[j2];
  2032.     x1i = -a[1] + a[j2 + 1];
  2033.     x2r = a[j1] + a[j3];
  2034.     x2i = a[j1 + 1] + a[j3 + 1];
  2035.     x3r = a[j1] - a[j3];
  2036.     x3i = a[j1 + 1] - a[j3 + 1];
  2037.     a[0] = x0r + x2r;
  2038.     a[1] = x0i - x2i;
  2039.     a[j1] = x0r - x2r;
  2040.     a[j1 + 1] = x0i + x2i;
  2041.     a[j2] = x1r + x3i;
  2042.     a[j2 + 1] = x1i + x3r;
  2043.     a[j3] = x1r - x3i;
  2044.     a[j3 + 1] = x1i - x3r;
  2045.     wn4r = w[1];
  2046.     csc1 = w[2];
  2047.     csc3 = w[3];
  2048.     wd1r = 1;
  2049.     wd1i = 0;
  2050.     wd3r = 1;
  2051.     wd3i = 0;
  2052.     k = 0;
  2053.     for (j = 2; j < mh - 2; j += 4) {
  2054.         k += 4;
  2055.         wk1r = csc1 * (wd1r + w[k]);
  2056.         wk1i = csc1 * (wd1i + w[k + 1]);
  2057.         wk3r = csc3 * (wd3r + w[k + 2]);
  2058.         wk3i = csc3 * (wd3i + w[k + 3]);
  2059.         wd1r = w[k];
  2060.         wd1i = w[k + 1];
  2061.         wd3r = w[k + 2];
  2062.         wd3i = w[k + 3];
  2063.         j1 = j + m;
  2064.         j2 = j1 + m;
  2065.         j3 = j2 + m;
  2066.         x0r = a[j] + a[j2];
  2067.         x0i = -a[j + 1] - a[j2 + 1];
  2068.         x1r = a[j] - a[j2];
  2069.         x1i = -a[j + 1] + a[j2 + 1];
  2070.         y0r = a[j + 2] + a[j2 + 2];
  2071.         y0i = -a[j + 3] - a[j2 + 3];
  2072.         y1r = a[j + 2] - a[j2 + 2];
  2073.         y1i = -a[j + 3] + a[j2 + 3];
  2074.         x2r = a[j1] + a[j3];
  2075.         x2i = a[j1 + 1] + a[j3 + 1];
  2076.         x3r = a[j1] - a[j3];
  2077.         x3i = a[j1 + 1] - a[j3 + 1];
  2078.         y2r = a[j1 + 2] + a[j3 + 2];
  2079.         y2i = a[j1 + 3] + a[j3 + 3];
  2080.         y3r = a[j1 + 2] - a[j3 + 2];
  2081.         y3i = a[j1 + 3] - a[j3 + 3];
  2082.         a[j] = x0r + x2r;
  2083.         a[j + 1] = x0i - x2i;
  2084.         a[j + 2] = y0r + y2r;
  2085.         a[j + 3] = y0i - y2i;
  2086.         a[j1] = x0r - x2r;
  2087.         a[j1 + 1] = x0i + x2i;
  2088.         a[j1 + 2] = y0r - y2r;
  2089.         a[j1 + 3] = y0i + y2i;
  2090.         x0r = x1r + x3i;
  2091.         x0i = x1i + x3r;
  2092.         a[j2] = wk1r * x0r - wk1i * x0i;
  2093.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  2094.         x0r = y1r + y3i;
  2095.         x0i = y1i + y3r;
  2096.         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
  2097.         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
  2098.         x0r = x1r - x3i;
  2099.         x0i = x1i - x3r;
  2100.         a[j3] = wk3r * x0r + wk3i * x0i;
  2101.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  2102.         x0r = y1r - y3i;
  2103.         x0i = y1i - y3r;
  2104.         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
  2105.         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
  2106.         j0 = m - j;
  2107.         j1 = j0 + m;
  2108.         j2 = j1 + m;
  2109.         j3 = j2 + m;
  2110.         x0r = a[j0] + a[j2];
  2111.         x0i = -a[j0 + 1] - a[j2 + 1];
  2112.         x1r = a[j0] - a[j2];
  2113.         x1i = -a[j0 + 1] + a[j2 + 1];
  2114.         y0r = a[j0 - 2] + a[j2 - 2];
  2115.         y0i = -a[j0 - 1] - a[j2 - 1];
  2116.         y1r = a[j0 - 2] - a[j2 - 2];
  2117.         y1i = -a[j0 - 1] + a[j2 - 1];
  2118.         x2r = a[j1] + a[j3];
  2119.         x2i = a[j1 + 1] + a[j3 + 1];
  2120.         x3r = a[j1] - a[j3];
  2121.         x3i = a[j1 + 1] - a[j3 + 1];
  2122.         y2r = a[j1 - 2] + a[j3 - 2];
  2123.         y2i = a[j1 - 1] + a[j3 - 1];
  2124.         y3r = a[j1 - 2] - a[j3 - 2];
  2125.         y3i = a[j1 - 1] - a[j3 - 1];
  2126.         a[j0] = x0r + x2r;
  2127.         a[j0 + 1] = x0i - x2i;
  2128.         a[j0 - 2] = y0r + y2r;
  2129.         a[j0 - 1] = y0i - y2i;
  2130.         a[j1] = x0r - x2r;
  2131.         a[j1 + 1] = x0i + x2i;
  2132.         a[j1 - 2] = y0r - y2r;
  2133.         a[j1 - 1] = y0i + y2i;
  2134.         x0r = x1r + x3i;
  2135.         x0i = x1i + x3r;
  2136.         a[j2] = wk1i * x0r - wk1r * x0i;
  2137.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  2138.         x0r = y1r + y3i;
  2139.         x0i = y1i + y3r;
  2140.         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
  2141.         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
  2142.         x0r = x1r - x3i;
  2143.         x0i = x1i - x3r;
  2144.         a[j3] = wk3i * x0r + wk3r * x0i;
  2145.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  2146.         x0r = y1r - y3i;
  2147.         x0i = y1i - y3r;
  2148.         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
  2149.         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
  2150.     }
  2151.     wk1r = csc1 * (wd1r + wn4r);
  2152.     wk1i = csc1 * (wd1i + wn4r);
  2153.     wk3r = csc3 * (wd3r - wn4r);
  2154.     wk3i = csc3 * (wd3i - wn4r);
  2155.     j0 = mh;
  2156.     j1 = j0 + m;
  2157.     j2 = j1 + m;
  2158.     j3 = j2 + m;
  2159.     x0r = a[j0 - 2] + a[j2 - 2];
  2160.     x0i = -a[j0 - 1] - a[j2 - 1];
  2161.     x1r = a[j0 - 2] - a[j2 - 2];
  2162.     x1i = -a[j0 - 1] + a[j2 - 1];
  2163.     x2r = a[j1 - 2] + a[j3 - 2];
  2164.     x2i = a[j1 - 1] + a[j3 - 1];
  2165.     x3r = a[j1 - 2] - a[j3 - 2];
  2166.     x3i = a[j1 - 1] - a[j3 - 1];
  2167.     a[j0 - 2] = x0r + x2r;
  2168.     a[j0 - 1] = x0i - x2i;
  2169.     a[j1 - 2] = x0r - x2r;
  2170.     a[j1 - 1] = x0i + x2i;
  2171.     x0r = x1r + x3i;
  2172.     x0i = x1i + x3r;
  2173.     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
  2174.     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
  2175.     x0r = x1r - x3i;
  2176.     x0i = x1i - x3r;
  2177.     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
  2178.     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
  2179.     x0r = a[j0] + a[j2];
  2180.     x0i = -a[j0 + 1] - a[j2 + 1];
  2181.     x1r = a[j0] - a[j2];
  2182.     x1i = -a[j0 + 1] + a[j2 + 1];
  2183.     x2r = a[j1] + a[j3];
  2184.     x2i = a[j1 + 1] + a[j3 + 1];
  2185.     x3r = a[j1] - a[j3];
  2186.     x3i = a[j1 + 1] - a[j3 + 1];
  2187.     a[j0] = x0r + x2r;
  2188.     a[j0 + 1] = x0i - x2i;
  2189.     a[j1] = x0r - x2r;
  2190.     a[j1 + 1] = x0i + x2i;
  2191.     x0r = x1r + x3i;
  2192.     x0i = x1i + x3r;
  2193.     a[j2] = wn4r * (x0r - x0i);
  2194.     a[j2 + 1] = wn4r * (x0i + x0r);
  2195.     x0r = x1r - x3i;
  2196.     x0i = x1i - x3r;
  2197.     a[j3] = -wn4r * (x0r + x0i);
  2198.     a[j3 + 1] = -wn4r * (x0i - x0r);
  2199.     x0r = a[j0 + 2] + a[j2 + 2];
  2200.     x0i = -a[j0 + 3] - a[j2 + 3];
  2201.     x1r = a[j0 + 2] - a[j2 + 2];
  2202.     x1i = -a[j0 + 3] + a[j2 + 3];
  2203.     x2r = a[j1 + 2] + a[j3 + 2];
  2204.     x2i = a[j1 + 3] + a[j3 + 3];
  2205.     x3r = a[j1 + 2] - a[j3 + 2];
  2206.     x3i = a[j1 + 3] - a[j3 + 3];
  2207.     a[j0 + 2] = x0r + x2r;
  2208.     a[j0 + 3] = x0i - x2i;
  2209.     a[j1 + 2] = x0r - x2r;
  2210.     a[j1 + 3] = x0i + x2i;
  2211.     x0r = x1r + x3i;
  2212.     x0i = x1i + x3r;
  2213.     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
  2214.     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
  2215.     x0r = x1r - x3i;
  2216.     x0i = x1i - x3r;
  2217.     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
  2218.     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
  2219. }
  2222. #ifdef USE_CDFT_THREADS
  2223. struct cdft_arg_st {
  2224.     int n0;
  2225.     int n;
  2226.     double *a;
  2227.     int nw;
  2228.     double *w;
  2229. };
  2230. typedef struct cdft_arg_st cdft_arg_t;
  2233. void cftrec4_th(int n, double *a, int nw, double *w)
  2234. {
  2235.     void *cftrec1_th(void *p);
  2236.     void *cftrec2_th(void *p);
  2237.     int i, idiv4, m, nthread;
  2238.     cdft_thread_t th[4];
  2239.     cdft_arg_t ag[4];
  2240.     
  2241.     nthread = 2;
  2242.     idiv4 = 0;
  2243.     m = n >> 1;
  2244.     if (n > CDFT_4THREADS_BEGIN_N) {
  2245.         nthread = 4;
  2246.         idiv4 = 1;
  2247.         m >>= 1;
  2248.     }
  2249.     for (i = 0; i < nthread; i++) {
  2250.         ag[i].n0 = n;
  2251.         ag[i].n = m;
  2252.         ag[i].a = &a[i * m];
  2253.         ag[i].nw = nw;
  2254.         ag[i].w = w;
  2255.         if (i != idiv4) {
  2256.             cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
  2257.         } else {
  2258.             cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
  2259.         }
  2260.     }
  2261.     for (i = 0; i < nthread; i++) {
  2262.         cdft_thread_wait(th[i]);
  2263.     }
  2264. }
  2267. void *cftrec1_th(void *p)
  2268. {
  2269.     int cfttree(int n, int j, int k, double *a, int nw, double *w);
  2270.     void cftleaf(int n, int isplt, double *a, int nw, double *w);
  2271.     void cftmdl1(int n, double *a, double *w);
  2272.     int isplt, j, k, m, n, n0, nw;
  2273.     double *a, *w;
  2274.     
  2275.     n0 = ((cdft_arg_t *) p)->n0;
  2276.     n = ((cdft_arg_t *) p)->n;
  2277.     a = ((cdft_arg_t *) p)->a;
  2278.     nw = ((cdft_arg_t *) p)->nw;
  2279.     w = ((cdft_arg_t *) p)->w;
  2280.     m = n0;
  2281.     while (m > 512) {
  2282.         m >>= 2;
  2283.         cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
  2284.     }
  2285.     cftleaf(m, 1, &a[n - m], nw, w);
  2286.     k = 0;
  2287.     for (j = n - m; j > 0; j -= m) {
  2288.         k++;
  2289.         isplt = cfttree(m, j, k, a, nw, w);
  2290.         cftleaf(m, isplt, &a[j - m], nw, w);
  2291.     }
  2292.     return (void *) 0;
  2293. }
  2296. void *cftrec2_th(void *p)
  2297. {
  2298.     int cfttree(int n, int j, int k, double *a, int nw, double *w);
  2299.     void cftleaf(int n, int isplt, double *a, int nw, double *w);
  2300.     void cftmdl2(int n, double *a, double *w);
  2301.     int isplt, j, k, m, n, n0, nw;
  2302.     double *a, *w;
  2303.     
  2304.     n0 = ((cdft_arg_t *) p)->n0;
  2305.     n = ((cdft_arg_t *) p)->n;
  2306.     a = ((cdft_arg_t *) p)->a;
  2307.     nw = ((cdft_arg_t *) p)->nw;
  2308.     w = ((cdft_arg_t *) p)->w;
  2309.     k = 1;
  2310.     m = n0;
  2311.     while (m > 512) {
  2312.         m >>= 2;
  2313.         k <<= 2;
  2314.         cftmdl2(m, &a[n - m], &w[nw - m]);
  2315.     }
  2316.     cftleaf(m, 0, &a[n - m], nw, w);
  2317.     k >>= 1;
  2318.     for (j = n - m; j > 0; j -= m) {
  2319.         k++;
  2320.         isplt = cfttree(m, j, k, a, nw, w);
  2321.         cftleaf(m, isplt, &a[j - m], nw, w);
  2322.     }
  2323.     return (void *) 0;
  2324. }
  2325. #endif /* USE_CDFT_THREADS */
  2328. void cftrec4(int n, double *a, int nw, double *w)
  2329. {
  2330.     int cfttree(int n, int j, int k, double *a, int nw, double *w);
  2331.     void cftleaf(int n, int isplt, double *a, int nw, double *w);
  2332.     void cftmdl1(int n, double *a, double *w);
  2333.     int isplt, j, k, m;
  2334.     
  2335.     m = n;
  2336.     while (m > 512) {
  2337.         m >>= 2;
  2338.         cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
  2339.     }
  2340.     cftleaf(m, 1, &a[n - m], nw, w);
  2341.     k = 0;
  2342.     for (j = n - m; j > 0; j -= m) {
  2343.         k++;
  2344.         isplt = cfttree(m, j, k, a, nw, w);
  2345.         cftleaf(m, isplt, &a[j - m], nw, w);
  2346.     }
  2347. }
  2350. int cfttree(int n, int j, int k, double *a, int nw, double *w)
  2351. {
  2352.     void cftmdl1(int n, double *a, double *w);
  2353.     void cftmdl2(int n, double *a, double *w);
  2354.     int i, isplt, m;
  2355.     
  2356.     if ((k & 3) != 0) {
  2357.         isplt = k & 1;
  2358.         if (isplt != 0) {
  2359.             cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
  2360.         } else {
  2361.             cftmdl2(n, &a[j - n], &w[nw - n]);
  2362.         }
  2363.     } else {
  2364.         m = n;
  2365.         for (i = k; (i & 3) == 0; i >>= 2) {
  2366.             m <<= 2;
  2367.         }
  2368.         isplt = i & 1;
  2369.         if (isplt != 0) {
  2370.             while (m > 128) {
  2371.                 cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
  2372.                 m >>= 2;
  2373.             }
  2374.         } else {
  2375.             while (m > 128) {
  2376.                 cftmdl2(m, &a[j - m], &w[nw - m]);
  2377.                 m >>= 2;
  2378.             }
  2379.         }
  2380.     }
  2381.     return isplt;
  2382. }
  2385. void cftleaf(int n, int isplt, double *a, int nw, double *w)
  2386. {
  2387.     void cftmdl1(int n, double *a, double *w);
  2388.     void cftmdl2(int n, double *a, double *w);
  2389.     void cftf161(double *a, double *w);
  2390.     void cftf162(double *a, double *w);
  2391.     void cftf081(double *a, double *w);
  2392.     void cftf082(double *a, double *w);
  2393.     
  2394.     if (n == 512) {
  2395.         cftmdl1(128, a, &w[nw - 64]);
  2396.         cftf161(a, &w[nw - 8]);
  2397.         cftf162(&a[32], &w[nw - 32]);
  2398.         cftf161(&a[64], &w[nw - 8]);
  2399.         cftf161(&a[96], &w[nw - 8]);
  2400.         cftmdl2(128, &a[128], &w[nw - 128]);
  2401.         cftf161(&a[128], &w[nw - 8]);
  2402.         cftf162(&a[160], &w[nw - 32]);
  2403.         cftf161(&a[192], &w[nw - 8]);
  2404.         cftf162(&a[224], &w[nw - 32]);
  2405.         cftmdl1(128, &a[256], &w[nw - 64]);
  2406.         cftf161(&a[256], &w[nw - 8]);
  2407.         cftf162(&a[288], &w[nw - 32]);
  2408.         cftf161(&a[320], &w[nw - 8]);
  2409.         cftf161(&a[352], &w[nw - 8]);
  2410.         if (isplt != 0) {
  2411.             cftmdl1(128, &a[384], &w[nw - 64]);
  2412.             cftf161(&a[480], &w[nw - 8]);
  2413.         } else {
  2414.             cftmdl2(128, &a[384], &w[nw - 128]);
  2415.             cftf162(&a[480], &w[nw - 32]);
  2416.         }
  2417.         cftf161(&a[384], &w[nw - 8]);
  2418.         cftf162(&a[416], &w[nw - 32]);
  2419.         cftf161(&a[448], &w[nw - 8]);
  2420.     } else {
  2421.         cftmdl1(64, a, &w[nw - 32]);
  2422.         cftf081(a, &w[nw - 8]);
  2423.         cftf082(&a[16], &w[nw - 8]);
  2424.         cftf081(&a[32], &w[nw - 8]);
  2425.         cftf081(&a[48], &w[nw - 8]);
  2426.         cftmdl2(64, &a[64], &w[nw - 64]);
  2427.         cftf081(&a[64], &w[nw - 8]);
  2428.         cftf082(&a[80], &w[nw - 8]);
  2429.         cftf081(&a[96], &w[nw - 8]);
  2430.         cftf082(&a[112], &w[nw - 8]);
  2431.         cftmdl1(64, &a[128], &w[nw - 32]);
  2432.         cftf081(&a[128], &w[nw - 8]);
  2433.         cftf082(&a[144], &w[nw - 8]);
  2434.         cftf081(&a[160], &w[nw - 8]);
  2435.         cftf081(&a[176], &w[nw - 8]);
  2436.         if (isplt != 0) {
  2437.             cftmdl1(64, &a[192], &w[nw - 32]);
  2438.             cftf081(&a[240], &w[nw - 8]);
  2439.         } else {
  2440.             cftmdl2(64, &a[192], &w[nw - 64]);
  2441.             cftf082(&a[240], &w[nw - 8]);
  2442.         }
  2443.         cftf081(&a[192], &w[nw - 8]);
  2444.         cftf082(&a[208], &w[nw - 8]);
  2445.         cftf081(&a[224], &w[nw - 8]);
  2446.     }
  2447. }
  2450. void cftmdl1(int n, double *a, double *w)
  2451. {
  2452.     int j, j0, j1, j2, j3, k, m, mh;
  2453.     double wn4r, wk1r, wk1i, wk3r, wk3i;
  2454.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  2455.     
  2456.     mh = n >> 3;
  2457.     m = 2 * mh;
  2458.     j1 = m;
  2459.     j2 = j1 + m;
  2460.     j3 = j2 + m;
  2461.     x0r = a[0] + a[j2];
  2462.     x0i = a[1] + a[j2 + 1];
  2463.     x1r = a[0] - a[j2];
  2464.     x1i = a[1] - a[j2 + 1];
  2465.     x2r = a[j1] + a[j3];
  2466.     x2i = a[j1 + 1] + a[j3 + 1];
  2467.     x3r = a[j1] - a[j3];
  2468.     x3i = a[j1 + 1] - a[j3 + 1];
  2469.     a[0] = x0r + x2r;
  2470.     a[1] = x0i + x2i;
  2471.     a[j1] = x0r - x2r;
  2472.     a[j1 + 1] = x0i - x2i;
  2473.     a[j2] = x1r - x3i;
  2474.     a[j2 + 1] = x1i + x3r;
  2475.     a[j3] = x1r + x3i;
  2476.     a[j3 + 1] = x1i - x3r;
  2477.     wn4r = w[1];
  2478.     k = 0;
  2479.     for (j = 2; j < mh; j += 2) {
  2480.         k += 4;
  2481.         wk1r = w[k];
  2482.         wk1i = w[k + 1];
  2483.         wk3r = w[k + 2];
  2484.         wk3i = w[k + 3];
  2485.         j1 = j + m;
  2486.         j2 = j1 + m;
  2487.         j3 = j2 + m;
  2488.         x0r = a[j] + a[j2];
  2489.         x0i = a[j + 1] + a[j2 + 1];
  2490.         x1r = a[j] - a[j2];
  2491.         x1i = a[j + 1] - a[j2 + 1];
  2492.         x2r = a[j1] + a[j3];
  2493.         x2i = a[j1 + 1] + a[j3 + 1];
  2494.         x3r = a[j1] - a[j3];
  2495.         x3i = a[j1 + 1] - a[j3 + 1];
  2496.         a[j] = x0r + x2r;
  2497.         a[j + 1] = x0i + x2i;
  2498.         a[j1] = x0r - x2r;
  2499.         a[j1 + 1] = x0i - x2i;
  2500.         x0r = x1r - x3i;
  2501.         x0i = x1i + x3r;
  2502.         a[j2] = wk1r * x0r - wk1i * x0i;
  2503.         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
  2504.         x0r = x1r + x3i;
  2505.         x0i = x1i - x3r;
  2506.         a[j3] = wk3r * x0r + wk3i * x0i;
  2507.         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
  2508.         j0 = m - j;
  2509.         j1 = j0 + m;
  2510.         j2 = j1 + m;
  2511.         j3 = j2 + m;
  2512.         x0r = a[j0] + a[j2];
  2513.         x0i = a[j0 + 1] + a[j2 + 1];
  2514.         x1r = a[j0] - a[j2];
  2515.         x1i = a[j0 + 1] - a[j2 + 1];
  2516.         x2r = a[j1] + a[j3];
  2517.         x2i = a[j1 + 1] + a[j3 + 1];
  2518.         x3r = a[j1] - a[j3];
  2519.         x3i = a[j1 + 1] - a[j3 + 1];
  2520.         a[j0] = x0r + x2r;
  2521.         a[j0 + 1] = x0i + x2i;
  2522.         a[j1] = x0r - x2r;
  2523.         a[j1 + 1] = x0i - x2i;
  2524.         x0r = x1r - x3i;
  2525.         x0i = x1i + x3r;
  2526.         a[j2] = wk1i * x0r - wk1r * x0i;
  2527.         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
  2528.         x0r = x1r + x3i;
  2529.         x0i = x1i - x3r;
  2530.         a[j3] = wk3i * x0r + wk3r * x0i;
  2531.         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
  2532.     }
  2533.     j0 = mh;
  2534.     j1 = j0 + m;
  2535.     j2 = j1 + m;
  2536.     j3 = j2 + m;
  2537.     x0r = a[j0] + a[j2];
  2538.     x0i = a[j0 + 1] + a[j2 + 1];
  2539.     x1r = a[j0] - a[j2];
  2540.     x1i = a[j0 + 1] - a[j2 + 1];
  2541.     x2r = a[j1] + a[j3];
  2542.     x2i = a[j1 + 1] + a[j3 + 1];
  2543.     x3r = a[j1] - a[j3];
  2544.     x3i = a[j1 + 1] - a[j3 + 1];
  2545.     a[j0] = x0r + x2r;
  2546.     a[j0 + 1] = x0i + x2i;
  2547.     a[j1] = x0r - x2r;
  2548.     a[j1 + 1] = x0i - x2i;
  2549.     x0r = x1r - x3i;
  2550.     x0i = x1i + x3r;
  2551.     a[j2] = wn4r * (x0r - x0i);
  2552.     a[j2 + 1] = wn4r * (x0i + x0r);
  2553.     x0r = x1r + x3i;
  2554.     x0i = x1i - x3r;
  2555.     a[j3] = -wn4r * (x0r + x0i);
  2556.     a[j3 + 1] = -wn4r * (x0i - x0r);
  2557. }
  2560. void cftmdl2(int n, double *a, double *w)
  2561. {
  2562.     int j, j0, j1, j2, j3, k, kr, m, mh;
  2563.     double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
  2564.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
  2565.     
  2566.     mh = n >> 3;
  2567.     m = 2 * mh;
  2568.     wn4r = w[1];
  2569.     j1 = m;
  2570.     j2 = j1 + m;
  2571.     j3 = j2 + m;
  2572.     x0r = a[0] - a[j2 + 1];
  2573.     x0i = a[1] + a[j2];
  2574.     x1r = a[0] + a[j2 + 1];
  2575.     x1i = a[1] - a[j2];
  2576.     x2r = a[j1] - a[j3 + 1];
  2577.     x2i = a[j1 + 1] + a[j3];
  2578.     x3r = a[j1] + a[j3 + 1];
  2579.     x3i = a[j1 + 1] - a[j3];
  2580.     y0r = wn4r * (x2r - x2i);
  2581.     y0i = wn4r * (x2i + x2r);
  2582.     a[0] = x0r + y0r;
  2583.     a[1] = x0i + y0i;
  2584.     a[j1] = x0r - y0r;
  2585.     a[j1 + 1] = x0i - y0i;
  2586.     y0r = wn4r * (x3r - x3i);
  2587.     y0i = wn4r * (x3i + x3r);
  2588.     a[j2] = x1r - y0i;
  2589.     a[j2 + 1] = x1i + y0r;
  2590.     a[j3] = x1r + y0i;
  2591.     a[j3 + 1] = x1i - y0r;
  2592.     k = 0;
  2593.     kr = 2 * m;
  2594.     for (j = 2; j < mh; j += 2) {
  2595.         k += 4;
  2596.         wk1r = w[k];
  2597.         wk1i = w[k + 1];
  2598.         wk3r = w[k + 2];
  2599.         wk3i = w[k + 3];
  2600.         kr -= 4;
  2601.         wd1i = w[kr];
  2602.         wd1r = w[kr + 1];
  2603.         wd3i = w[kr + 2];
  2604.         wd3r = w[kr + 3];
  2605.         j1 = j + m;
  2606.         j2 = j1 + m;
  2607.         j3 = j2 + m;
  2608.         x0r = a[j] - a[j2 + 1];
  2609.         x0i = a[j + 1] + a[j2];
  2610.         x1r = a[j] + a[j2 + 1];
  2611.         x1i = a[j + 1] - a[j2];
  2612.         x2r = a[j1] - a[j3 + 1];
  2613.         x2i = a[j1 + 1] + a[j3];
  2614.         x3r = a[j1] + a[j3 + 1];
  2615.         x3i = a[j1 + 1] - a[j3];
  2616.         y0r = wk1r * x0r - wk1i * x0i;
  2617.         y0i = wk1r * x0i + wk1i * x0r;
  2618.         y2r = wd1r * x2r - wd1i * x2i;
  2619.         y2i = wd1r * x2i + wd1i * x2r;
  2620.         a[j] = y0r + y2r;
  2621.         a[j + 1] = y0i + y2i;
  2622.         a[j1] = y0r - y2r;
  2623.         a[j1 + 1] = y0i - y2i;
  2624.         y0r = wk3r * x1r + wk3i * x1i;
  2625.         y0i = wk3r * x1i - wk3i * x1r;
  2626.         y2r = wd3r * x3r + wd3i * x3i;
  2627.         y2i = wd3r * x3i - wd3i * x3r;
  2628.         a[j2] = y0r + y2r;
  2629.         a[j2 + 1] = y0i + y2i;
  2630.         a[j3] = y0r - y2r;
  2631.         a[j3 + 1] = y0i - y2i;
  2632.         j0 = m - j;
  2633.         j1 = j0 + m;
  2634.         j2 = j1 + m;
  2635.         j3 = j2 + m;
  2636.         x0r = a[j0] - a[j2 + 1];
  2637.         x0i = a[j0 + 1] + a[j2];
  2638.         x1r = a[j0] + a[j2 + 1];
  2639.         x1i = a[j0 + 1] - a[j2];
  2640.         x2r = a[j1] - a[j3 + 1];
  2641.         x2i = a[j1 + 1] + a[j3];
  2642.         x3r = a[j1] + a[j3 + 1];
  2643.         x3i = a[j1 + 1] - a[j3];
  2644.         y0r = wd1i * x0r - wd1r * x0i;
  2645.         y0i = wd1i * x0i + wd1r * x0r;
  2646.         y2r = wk1i * x2r - wk1r * x2i;
  2647.         y2i = wk1i * x2i + wk1r * x2r;
  2648.         a[j0] = y0r + y2r;
  2649.         a[j0 + 1] = y0i + y2i;
  2650.         a[j1] = y0r - y2r;
  2651.         a[j1 + 1] = y0i - y2i;
  2652.         y0r = wd3i * x1r + wd3r * x1i;
  2653.         y0i = wd3i * x1i - wd3r * x1r;
  2654.         y2r = wk3i * x3r + wk3r * x3i;
  2655.         y2i = wk3i * x3i - wk3r * x3r;
  2656.         a[j2] = y0r + y2r;
  2657.         a[j2 + 1] = y0i + y2i;
  2658.         a[j3] = y0r - y2r;
  2659.         a[j3 + 1] = y0i - y2i;
  2660.     }
  2661.     wk1r = w[m];
  2662.     wk1i = w[m + 1];
  2663.     j0 = mh;
  2664.     j1 = j0 + m;
  2665.     j2 = j1 + m;
  2666.     j3 = j2 + m;
  2667.     x0r = a[j0] - a[j2 + 1];
  2668.     x0i = a[j0 + 1] + a[j2];
  2669.     x1r = a[j0] + a[j2 + 1];
  2670.     x1i = a[j0 + 1] - a[j2];
  2671.     x2r = a[j1] - a[j3 + 1];
  2672.     x2i = a[j1 + 1] + a[j3];
  2673.     x3r = a[j1] + a[j3 + 1];
  2674.     x3i = a[j1 + 1] - a[j3];
  2675.     y0r = wk1r * x0r - wk1i * x0i;
  2676.     y0i = wk1r * x0i + wk1i * x0r;
  2677.     y2r = wk1i * x2r - wk1r * x2i;
  2678.     y2i = wk1i * x2i + wk1r * x2r;
  2679.     a[j0] = y0r + y2r;
  2680.     a[j0 + 1] = y0i + y2i;
  2681.     a[j1] = y0r - y2r;
  2682.     a[j1 + 1] = y0i - y2i;
  2683.     y0r = wk1i * x1r - wk1r * x1i;
  2684.     y0i = wk1i * x1i + wk1r * x1r;
  2685.     y2r = wk1r * x3r - wk1i * x3i;
  2686.     y2i = wk1r * x3i + wk1i * x3r;
  2687.     a[j2] = y0r - y2r;
  2688.     a[j2 + 1] = y0i - y2i;
  2689.     a[j3] = y0r + y2r;
  2690.     a[j3 + 1] = y0i + y2i;
  2691. }
  2694. void cftfx41(int n, double *a, int nw, double *w)
  2695. {
  2696.     void cftf161(double *a, double *w);
  2697.     void cftf162(double *a, double *w);
  2698.     void cftf081(double *a, double *w);
  2699.     void cftf082(double *a, double *w);
  2700.     
  2701.     if (n == 128) {
  2702.         cftf161(a, &w[nw - 8]);
  2703.         cftf162(&a[32], &w[nw - 32]);
  2704.         cftf161(&a[64], &w[nw - 8]);
  2705.         cftf161(&a[96], &w[nw - 8]);
  2706.     } else {
  2707.         cftf081(a, &w[nw - 8]);
  2708.         cftf082(&a[16], &w[nw - 8]);
  2709.         cftf081(&a[32], &w[nw - 8]);
  2710.         cftf081(&a[48], &w[nw - 8]);
  2711.     }
  2712. }
  2715. void cftf161(double *a, double *w)
  2716. {
  2717.     double wn4r, wk1r, wk1i, 
  2718.         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  2719.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2720.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
  2721.         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
  2722.         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
  2723.     
  2724.     wn4r = w[1];
  2725.     wk1r = w[2];
  2726.     wk1i = w[3];
  2727.     x0r = a[0] + a[16];
  2728.     x0i = a[1] + a[17];
  2729.     x1r = a[0] - a[16];
  2730.     x1i = a[1] - a[17];
  2731.     x2r = a[8] + a[24];
  2732.     x2i = a[9] + a[25];
  2733.     x3r = a[8] - a[24];
  2734.     x3i = a[9] - a[25];
  2735.     y0r = x0r + x2r;
  2736.     y0i = x0i + x2i;
  2737.     y4r = x0r - x2r;
  2738.     y4i = x0i - x2i;
  2739.     y8r = x1r - x3i;
  2740.     y8i = x1i + x3r;
  2741.     y12r = x1r + x3i;
  2742.     y12i = x1i - x3r;
  2743.     x0r = a[2] + a[18];
  2744.     x0i = a[3] + a[19];
  2745.     x1r = a[2] - a[18];
  2746.     x1i = a[3] - a[19];
  2747.     x2r = a[10] + a[26];
  2748.     x2i = a[11] + a[27];
  2749.     x3r = a[10] - a[26];
  2750.     x3i = a[11] - a[27];
  2751.     y1r = x0r + x2r;
  2752.     y1i = x0i + x2i;
  2753.     y5r = x0r - x2r;
  2754.     y5i = x0i - x2i;
  2755.     x0r = x1r - x3i;
  2756.     x0i = x1i + x3r;
  2757.     y9r = wk1r * x0r - wk1i * x0i;
  2758.     y9i = wk1r * x0i + wk1i * x0r;
  2759.     x0r = x1r + x3i;
  2760.     x0i = x1i - x3r;
  2761.     y13r = wk1i * x0r - wk1r * x0i;
  2762.     y13i = wk1i * x0i + wk1r * x0r;
  2763.     x0r = a[4] + a[20];
  2764.     x0i = a[5] + a[21];
  2765.     x1r = a[4] - a[20];
  2766.     x1i = a[5] - a[21];
  2767.     x2r = a[12] + a[28];
  2768.     x2i = a[13] + a[29];
  2769.     x3r = a[12] - a[28];
  2770.     x3i = a[13] - a[29];
  2771.     y2r = x0r + x2r;
  2772.     y2i = x0i + x2i;
  2773.     y6r = x0r - x2r;
  2774.     y6i = x0i - x2i;
  2775.     x0r = x1r - x3i;
  2776.     x0i = x1i + x3r;
  2777.     y10r = wn4r * (x0r - x0i);
  2778.     y10i = wn4r * (x0i + x0r);
  2779.     x0r = x1r + x3i;
  2780.     x0i = x1i - x3r;
  2781.     y14r = wn4r * (x0r + x0i);
  2782.     y14i = wn4r * (x0i - x0r);
  2783.     x0r = a[6] + a[22];
  2784.     x0i = a[7] + a[23];
  2785.     x1r = a[6] - a[22];
  2786.     x1i = a[7] - a[23];
  2787.     x2r = a[14] + a[30];
  2788.     x2i = a[15] + a[31];
  2789.     x3r = a[14] - a[30];
  2790.     x3i = a[15] - a[31];
  2791.     y3r = x0r + x2r;
  2792.     y3i = x0i + x2i;
  2793.     y7r = x0r - x2r;
  2794.     y7i = x0i - x2i;
  2795.     x0r = x1r - x3i;
  2796.     x0i = x1i + x3r;
  2797.     y11r = wk1i * x0r - wk1r * x0i;
  2798.     y11i = wk1i * x0i + wk1r * x0r;
  2799.     x0r = x1r + x3i;
  2800.     x0i = x1i - x3r;
  2801.     y15r = wk1r * x0r - wk1i * x0i;
  2802.     y15i = wk1r * x0i + wk1i * x0r;
  2803.     x0r = y12r - y14r;
  2804.     x0i = y12i - y14i;
  2805.     x1r = y12r + y14r;
  2806.     x1i = y12i + y14i;
  2807.     x2r = y13r - y15r;
  2808.     x2i = y13i - y15i;
  2809.     x3r = y13r + y15r;
  2810.     x3i = y13i + y15i;
  2811.     a[24] = x0r + x2r;
  2812.     a[25] = x0i + x2i;
  2813.     a[26] = x0r - x2r;
  2814.     a[27] = x0i - x2i;
  2815.     a[28] = x1r - x3i;
  2816.     a[29] = x1i + x3r;
  2817.     a[30] = x1r + x3i;
  2818.     a[31] = x1i - x3r;
  2819.     x0r = y8r + y10r;
  2820.     x0i = y8i + y10i;
  2821.     x1r = y8r - y10r;
  2822.     x1i = y8i - y10i;
  2823.     x2r = y9r + y11r;
  2824.     x2i = y9i + y11i;
  2825.     x3r = y9r - y11r;
  2826.     x3i = y9i - y11i;
  2827.     a[16] = x0r + x2r;
  2828.     a[17] = x0i + x2i;
  2829.     a[18] = x0r - x2r;
  2830.     a[19] = x0i - x2i;
  2831.     a[20] = x1r - x3i;
  2832.     a[21] = x1i + x3r;
  2833.     a[22] = x1r + x3i;
  2834.     a[23] = x1i - x3r;
  2835.     x0r = y5r - y7i;
  2836.     x0i = y5i + y7r;
  2837.     x2r = wn4r * (x0r - x0i);
  2838.     x2i = wn4r * (x0i + x0r);
  2839.     x0r = y5r + y7i;
  2840.     x0i = y5i - y7r;
  2841.     x3r = wn4r * (x0r - x0i);
  2842.     x3i = wn4r * (x0i + x0r);
  2843.     x0r = y4r - y6i;
  2844.     x0i = y4i + y6r;
  2845.     x1r = y4r + y6i;
  2846.     x1i = y4i - y6r;
  2847.     a[8] = x0r + x2r;
  2848.     a[9] = x0i + x2i;
  2849.     a[10] = x0r - x2r;
  2850.     a[11] = x0i - x2i;
  2851.     a[12] = x1r - x3i;
  2852.     a[13] = x1i + x3r;
  2853.     a[14] = x1r + x3i;
  2854.     a[15] = x1i - x3r;
  2855.     x0r = y0r + y2r;
  2856.     x0i = y0i + y2i;
  2857.     x1r = y0r - y2r;
  2858.     x1i = y0i - y2i;
  2859.     x2r = y1r + y3r;
  2860.     x2i = y1i + y3i;
  2861.     x3r = y1r - y3r;
  2862.     x3i = y1i - y3i;
  2863.     a[0] = x0r + x2r;
  2864.     a[1] = x0i + x2i;
  2865.     a[2] = x0r - x2r;
  2866.     a[3] = x0i - x2i;
  2867.     a[4] = x1r - x3i;
  2868.     a[5] = x1i + x3r;
  2869.     a[6] = x1r + x3i;
  2870.     a[7] = x1i - x3r;
  2871. }
  2874. void cftf162(double *a, double *w)
  2875. {
  2876.     double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
  2877.         x0r, x0i, x1r, x1i, x2r, x2i, 
  2878.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  2879.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i, 
  2880.         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i, 
  2881.         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
  2882.     
  2883.     wn4r = w[1];
  2884.     wk1r = w[4];
  2885.     wk1i = w[5];
  2886.     wk3r = w[6];
  2887.     wk3i = -w[7];
  2888.     wk2r = w[8];
  2889.     wk2i = w[9];
  2890.     x1r = a[0] - a[17];
  2891.     x1i = a[1] + a[16];
  2892.     x0r = a[8] - a[25];
  2893.     x0i = a[9] + a[24];
  2894.     x2r = wn4r * (x0r - x0i);
  2895.     x2i = wn4r * (x0i + x0r);
  2896.     y0r = x1r + x2r;
  2897.     y0i = x1i + x2i;
  2898.     y4r = x1r - x2r;
  2899.     y4i = x1i - x2i;
  2900.     x1r = a[0] + a[17];
  2901.     x1i = a[1] - a[16];
  2902.     x0r = a[8] + a[25];
  2903.     x0i = a[9] - a[24];
  2904.     x2r = wn4r * (x0r - x0i);
  2905.     x2i = wn4r * (x0i + x0r);
  2906.     y8r = x1r - x2i;
  2907.     y8i = x1i + x2r;
  2908.     y12r = x1r + x2i;
  2909.     y12i = x1i - x2r;
  2910.     x0r = a[2] - a[19];
  2911.     x0i = a[3] + a[18];
  2912.     x1r = wk1r * x0r - wk1i * x0i;
  2913.     x1i = wk1r * x0i + wk1i * x0r;
  2914.     x0r = a[10] - a[27];
  2915.     x0i = a[11] + a[26];
  2916.     x2r = wk3i * x0r - wk3r * x0i;
  2917.     x2i = wk3i * x0i + wk3r * x0r;
  2918.     y1r = x1r + x2r;
  2919.     y1i = x1i + x2i;
  2920.     y5r = x1r - x2r;
  2921.     y5i = x1i - x2i;
  2922.     x0r = a[2] + a[19];
  2923.     x0i = a[3] - a[18];
  2924.     x1r = wk3r * x0r - wk3i * x0i;
  2925.     x1i = wk3r * x0i + wk3i * x0r;
  2926.     x0r = a[10] + a[27];
  2927.     x0i = a[11] - a[26];
  2928.     x2r = wk1r * x0r + wk1i * x0i;
  2929.     x2i = wk1r * x0i - wk1i * x0r;
  2930.     y9r = x1r - x2r;
  2931.     y9i = x1i - x2i;
  2932.     y13r = x1r + x2r;
  2933.     y13i = x1i + x2i;
  2934.     x0r = a[4] - a[21];
  2935.     x0i = a[5] + a[20];
  2936.     x1r = wk2r * x0r - wk2i * x0i;
  2937.     x1i = wk2r * x0i + wk2i * x0r;
  2938.     x0r = a[12] - a[29];
  2939.     x0i = a[13] + a[28];
  2940.     x2r = wk2i * x0r - wk2r * x0i;
  2941.     x2i = wk2i * x0i + wk2r * x0r;
  2942.     y2r = x1r + x2r;
  2943.     y2i = x1i + x2i;
  2944.     y6r = x1r - x2r;
  2945.     y6i = x1i - x2i;
  2946.     x0r = a[4] + a[21];
  2947.     x0i = a[5] - a[20];
  2948.     x1r = wk2i * x0r - wk2r * x0i;
  2949.     x1i = wk2i * x0i + wk2r * x0r;
  2950.     x0r = a[12] + a[29];
  2951.     x0i = a[13] - a[28];
  2952.     x2r = wk2r * x0r - wk2i * x0i;
  2953.     x2i = wk2r * x0i + wk2i * x0r;
  2954.     y10r = x1r - x2r;
  2955.     y10i = x1i - x2i;
  2956.     y14r = x1r + x2r;
  2957.     y14i = x1i + x2i;
  2958.     x0r = a[6] - a[23];
  2959.     x0i = a[7] + a[22];
  2960.     x1r = wk3r * x0r - wk3i * x0i;
  2961.     x1i = wk3r * x0i + wk3i * x0r;
  2962.     x0r = a[14] - a[31];
  2963.     x0i = a[15] + a[30];
  2964.     x2r = wk1i * x0r - wk1r * x0i;
  2965.     x2i = wk1i * x0i + wk1r * x0r;
  2966.     y3r = x1r + x2r;
  2967.     y3i = x1i + x2i;
  2968.     y7r = x1r - x2r;
  2969.     y7i = x1i - x2i;
  2970.     x0r = a[6] + a[23];
  2971.     x0i = a[7] - a[22];
  2972.     x1r = wk1i * x0r + wk1r * x0i;
  2973.     x1i = wk1i * x0i - wk1r * x0r;
  2974.     x0r = a[14] + a[31];
  2975.     x0i = a[15] - a[30];
  2976.     x2r = wk3i * x0r - wk3r * x0i;
  2977.     x2i = wk3i * x0i + wk3r * x0r;
  2978.     y11r = x1r + x2r;
  2979.     y11i = x1i + x2i;
  2980.     y15r = x1r - x2r;
  2981.     y15i = x1i - x2i;
  2982.     x1r = y0r + y2r;
  2983.     x1i = y0i + y2i;
  2984.     x2r = y1r + y3r;
  2985.     x2i = y1i + y3i;
  2986.     a[0] = x1r + x2r;
  2987.     a[1] = x1i + x2i;
  2988.     a[2] = x1r - x2r;
  2989.     a[3] = x1i - x2i;
  2990.     x1r = y0r - y2r;
  2991.     x1i = y0i - y2i;
  2992.     x2r = y1r - y3r;
  2993.     x2i = y1i - y3i;
  2994.     a[4] = x1r - x2i;
  2995.     a[5] = x1i + x2r;
  2996.     a[6] = x1r + x2i;
  2997.     a[7] = x1i - x2r;
  2998.     x1r = y4r - y6i;
  2999.     x1i = y4i + y6r;
  3000.     x0r = y5r - y7i;
  3001.     x0i = y5i + y7r;
  3002.     x2r = wn4r * (x0r - x0i);
  3003.     x2i = wn4r * (x0i + x0r);
  3004.     a[8] = x1r + x2r;
  3005.     a[9] = x1i + x2i;
  3006.     a[10] = x1r - x2r;
  3007.     a[11] = x1i - x2i;
  3008.     x1r = y4r + y6i;
  3009.     x1i = y4i - y6r;
  3010.     x0r = y5r + y7i;
  3011.     x0i = y5i - y7r;
  3012.     x2r = wn4r * (x0r - x0i);
  3013.     x2i = wn4r * (x0i + x0r);
  3014.     a[12] = x1r - x2i;
  3015.     a[13] = x1i + x2r;
  3016.     a[14] = x1r + x2i;
  3017.     a[15] = x1i - x2r;
  3018.     x1r = y8r + y10r;
  3019.     x1i = y8i + y10i;
  3020.     x2r = y9r - y11r;
  3021.     x2i = y9i - y11i;
  3022.     a[16] = x1r + x2r;
  3023.     a[17] = x1i + x2i;
  3024.     a[18] = x1r - x2r;
  3025.     a[19] = x1i - x2i;
  3026.     x1r = y8r - y10r;
  3027.     x1i = y8i - y10i;
  3028.     x2r = y9r + y11r;
  3029.     x2i = y9i + y11i;
  3030.     a[20] = x1r - x2i;
  3031.     a[21] = x1i + x2r;
  3032.     a[22] = x1r + x2i;
  3033.     a[23] = x1i - x2r;
  3034.     x1r = y12r - y14i;
  3035.     x1i = y12i + y14r;
  3036.     x0r = y13r + y15i;
  3037.     x0i = y13i - y15r;
  3038.     x2r = wn4r * (x0r - x0i);
  3039.     x2i = wn4r * (x0i + x0r);
  3040.     a[24] = x1r + x2r;
  3041.     a[25] = x1i + x2i;
  3042.     a[26] = x1r - x2r;
  3043.     a[27] = x1i - x2i;
  3044.     x1r = y12r + y14i;
  3045.     x1i = y12i - y14r;
  3046.     x0r = y13r - y15i;
  3047.     x0i = y13i + y15r;
  3048.     x2r = wn4r * (x0r - x0i);
  3049.     x2i = wn4r * (x0i + x0r);
  3050.     a[28] = x1r - x2i;
  3051.     a[29] = x1i + x2r;
  3052.     a[30] = x1r + x2i;
  3053.     a[31] = x1i - x2r;
  3054. }
  3057. void cftf081(double *a, double *w)
  3058. {
  3059.     double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
  3060.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  3061.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
  3062.     
  3063.     wn4r = w[1];
  3064.     x0r = a[0] + a[8];
  3065.     x0i = a[1] + a[9];
  3066.     x1r = a[0] - a[8];
  3067.     x1i = a[1] - a[9];
  3068.     x2r = a[4] + a[12];
  3069.     x2i = a[5] + a[13];
  3070.     x3r = a[4] - a[12];
  3071.     x3i = a[5] - a[13];
  3072.     y0r = x0r + x2r;
  3073.     y0i = x0i + x2i;
  3074.     y2r = x0r - x2r;
  3075.     y2i = x0i - x2i;
  3076.     y1r = x1r - x3i;
  3077.     y1i = x1i + x3r;
  3078.     y3r = x1r + x3i;
  3079.     y3i = x1i - x3r;
  3080.     x0r = a[2] + a[10];
  3081.     x0i = a[3] + a[11];
  3082.     x1r = a[2] - a[10];
  3083.     x1i = a[3] - a[11];
  3084.     x2r = a[6] + a[14];
  3085.     x2i = a[7] + a[15];
  3086.     x3r = a[6] - a[14];
  3087.     x3i = a[7] - a[15];
  3088.     y4r = x0r + x2r;
  3089.     y4i = x0i + x2i;
  3090.     y6r = x0r - x2r;
  3091.     y6i = x0i - x2i;
  3092.     x0r = x1r - x3i;
  3093.     x0i = x1i + x3r;
  3094.     x2r = x1r + x3i;
  3095.     x2i = x1i - x3r;
  3096.     y5r = wn4r * (x0r - x0i);
  3097.     y5i = wn4r * (x0r + x0i);
  3098.     y7r = wn4r * (x2r - x2i);
  3099.     y7i = wn4r * (x2r + x2i);
  3100.     a[8] = y1r + y5r;
  3101.     a[9] = y1i + y5i;
  3102.     a[10] = y1r - y5r;
  3103.     a[11] = y1i - y5i;
  3104.     a[12] = y3r - y7i;
  3105.     a[13] = y3i + y7r;
  3106.     a[14] = y3r + y7i;
  3107.     a[15] = y3i - y7r;
  3108.     a[0] = y0r + y4r;
  3109.     a[1] = y0i + y4i;
  3110.     a[2] = y0r - y4r;
  3111.     a[3] = y0i - y4i;
  3112.     a[4] = y2r - y6i;
  3113.     a[5] = y2i + y6r;
  3114.     a[6] = y2r + y6i;
  3115.     a[7] = y2i - y6r;
  3116. }
  3119. void cftf082(double *a, double *w)
  3120. {
  3121.     double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
  3122.         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
  3123.         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
  3124.     
  3125.     wn4r = w[1];
  3126.     wk1r = w[2];
  3127.     wk1i = w[3];
  3128.     y0r = a[0] - a[9];
  3129.     y0i = a[1] + a[8];
  3130.     y1r = a[0] + a[9];
  3131.     y1i = a[1] - a[8];
  3132.     x0r = a[4] - a[13];
  3133.     x0i = a[5] + a[12];
  3134.     y2r = wn4r * (x0r - x0i);
  3135.     y2i = wn4r * (x0i + x0r);
  3136.     x0r = a[4] + a[13];
  3137.     x0i = a[5] - a[12];
  3138.     y3r = wn4r * (x0r - x0i);
  3139.     y3i = wn4r * (x0i + x0r);
  3140.     x0r = a[2] - a[11];
  3141.     x0i = a[3] + a[10];
  3142.     y4r = wk1r * x0r - wk1i * x0i;
  3143.     y4i = wk1r * x0i + wk1i * x0r;
  3144.     x0r = a[2] + a[11];
  3145.     x0i = a[3] - a[10];
  3146.     y5r = wk1i * x0r - wk1r * x0i;
  3147.     y5i = wk1i * x0i + wk1r * x0r;
  3148.     x0r = a[6] - a[15];
  3149.     x0i = a[7] + a[14];
  3150.     y6r = wk1i * x0r - wk1r * x0i;
  3151.     y6i = wk1i * x0i + wk1r * x0r;
  3152.     x0r = a[6] + a[15];
  3153.     x0i = a[7] - a[14];
  3154.     y7r = wk1r * x0r - wk1i * x0i;
  3155.     y7i = wk1r * x0i + wk1i * x0r;
  3156.     x0r = y0r + y2r;
  3157.     x0i = y0i + y2i;
  3158.     x1r = y4r + y6r;
  3159.     x1i = y4i + y6i;
  3160.     a[0] = x0r + x1r;
  3161.     a[1] = x0i + x1i;
  3162.     a[2] = x0r - x1r;
  3163.     a[3] = x0i - x1i;
  3164.     x0r = y0r - y2r;
  3165.     x0i = y0i - y2i;
  3166.     x1r = y4r - y6r;
  3167.     x1i = y4i - y6i;
  3168.     a[4] = x0r - x1i;
  3169.     a[5] = x0i + x1r;
  3170.     a[6] = x0r + x1i;
  3171.     a[7] = x0i - x1r;
  3172.     x0r = y1r - y3i;
  3173.     x0i = y1i + y3r;
  3174.     x1r = y5r - y7r;
  3175.     x1i = y5i - y7i;
  3176.     a[8] = x0r + x1r;
  3177.     a[9] = x0i + x1i;
  3178.     a[10] = x0r - x1r;
  3179.     a[11] = x0i - x1i;
  3180.     x0r = y1r + y3i;
  3181.     x0i = y1i - y3r;
  3182.     x1r = y5r + y7r;
  3183.     x1i = y5i + y7i;
  3184.     a[12] = x0r - x1i;
  3185.     a[13] = x0i + x1r;
  3186.     a[14] = x0r + x1i;
  3187.     a[15] = x0i - x1r;
  3188. }
  3191. void cftf040(double *a)
  3192. {
  3193.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  3194.     
  3195.     x0r = a[0] + a[4];
  3196.     x0i = a[1] + a[5];
  3197.     x1r = a[0] - a[4];
  3198.     x1i = a[1] - a[5];
  3199.     x2r = a[2] + a[6];
  3200.     x2i = a[3] + a[7];
  3201.     x3r = a[2] - a[6];
  3202.     x3i = a[3] - a[7];
  3203.     a[0] = x0r + x2r;
  3204.     a[1] = x0i + x2i;
  3205.     a[2] = x1r - x3i;
  3206.     a[3] = x1i + x3r;
  3207.     a[4] = x0r - x2r;
  3208.     a[5] = x0i - x2i;
  3209.     a[6] = x1r + x3i;
  3210.     a[7] = x1i - x3r;
  3211. }
  3214. void cftb040(double *a)
  3215. {
  3216.     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  3217.     
  3218.     x0r = a[0] + a[4];
  3219.     x0i = a[1] + a[5];
  3220.     x1r = a[0] - a[4];
  3221.     x1i = a[1] - a[5];
  3222.     x2r = a[2] + a[6];
  3223.     x2i = a[3] + a[7];
  3224.     x3r = a[2] - a[6];
  3225.     x3i = a[3] - a[7];
  3226.     a[0] = x0r + x2r;
  3227.     a[1] = x0i + x2i;
  3228.     a[2] = x1r + x3i;
  3229.     a[3] = x1i - x3r;
  3230.     a[4] = x0r - x2r;
  3231.     a[5] = x0i - x2i;
  3232.     a[6] = x1r - x3i;
  3233.     a[7] = x1i + x3r;
  3234. }
  3237. void cftx020(double *a)
  3238. {
  3239.     double x0r, x0i;
  3240.     
  3241.     x0r = a[0] - a[2];
  3242.     x0i = a[1] - a[3];
  3243.     a[0] += a[2];
  3244.     a[1] += a[3];
  3245.     a[2] = x0r;
  3246.     a[3] = x0i;
  3247. }
  3250. void rftfsub(int n, double *a, int nc, double *c)
  3251. {
  3252.     int j, k, kk, ks, m;
  3253.     double wkr, wki, xr, xi, yr, yi;
  3254.     
  3255.     m = n >> 1;
  3256.     ks = 2 * nc / m;
  3257.     kk = 0;
  3258.     for (j = 2; j < m; j += 2) {
  3259.         k = n - j;
  3260.         kk += ks;
  3261.         wkr = 0.5 - c[nc - kk];
  3262.         wki = c[kk];
  3263.         xr = a[j] - a[k];
  3264.         xi = a[j + 1] + a[k + 1];
  3265.         yr = wkr * xr - wki * xi;
  3266.         yi = wkr * xi + wki * xr;
  3267.         a[j] -= yr;
  3268.         a[j + 1] -= yi;
  3269.         a[k] += yr;
  3270.         a[k + 1] -= yi;
  3271.     }
  3272. }
  3275. void rftbsub(int n, double *a, int nc, double *c)
  3276. {
  3277.     int j, k, kk, ks, m;
  3278.     double wkr, wki, xr, xi, yr, yi;
  3279.     
  3280.     m = n >> 1;
  3281.     ks = 2 * nc / m;
  3282.     kk = 0;
  3283.     for (j = 2; j < m; j += 2) {
  3284.         k = n - j;
  3285.         kk += ks;
  3286.         wkr = 0.5 - c[nc - kk];
  3287.         wki = c[kk];
  3288.         xr = a[j] - a[k];
  3289.         xi = a[j + 1] + a[k + 1];
  3290.         yr = wkr * xr + wki * xi;
  3291.         yi = wkr * xi - wki * xr;
  3292.         a[j] -= yr;
  3293.         a[j + 1] -= yi;
  3294.         a[k] += yr;
  3295.         a[k + 1] -= yi;
  3296.     }
  3297. }
  3300. void dctsub(int n, double *a, int nc, double *c)
  3301. {
  3302.     int j, k, kk, ks, m;
  3303.     double wkr, wki, xr;
  3304.     
  3305.     m = n >> 1;
  3306.     ks = nc / n;
  3307.     kk = 0;
  3308.     for (j = 1; j < m; j++) {
  3309.         k = n - j;
  3310.         kk += ks;
  3311.         wkr = c[kk] - c[nc - kk];
  3312.         wki = c[kk] + c[nc - kk];
  3313.         xr = wki * a[j] - wkr * a[k];
  3314.         a[j] = wkr * a[j] + wki * a[k];
  3315.         a[k] = xr;
  3316.     }
  3317.     a[m] *= c[0];
  3318. }
  3321. void dstsub(int n, double *a, int nc, double *c)
  3322. {
  3323.     int j, k, kk, ks, m;
  3324.     double wkr, wki, xr;
  3325.     
  3326.     m = n >> 1;
  3327.     ks = nc / n;
  3328.     kk = 0;
  3329.     for (j = 1; j < m; j++) {
  3330.         k = n - j;
  3331.         kk += ks;
  3332.         wkr = c[kk] - c[nc - kk];
  3333.         wki = c[kk] + c[nc - kk];
  3334.         xr = wki * a[k] - wkr * a[j];
  3335.         a[k] = wkr * a[k] + wki * a[j];
  3336.         a[j] = xr;
  3337.     }
  3338.     a[m] *= c[0];
  3339. }