git.sesse.net Git - vlc/blob - modules/visualization/galaktos/fftsg.c

   1 /*****************************************************************************
   2  * fftsg.c:
   3  *****************************************************************************
   4  * Copyright (C) 2004 the VideoLAN team
   5  * $Id$
   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  *****************************************************************************/
  24
  25
  26
  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
  56
  57
  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         .
  99
 100
 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         .
 148
 149
 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         .
 187
 188
 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         .
 234
 235
 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         .
 270
 271
 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         .
 303
 304
 305 Appendix :
 306     The cos/sin table is recalculated when the larger table required.
 307     w[] and ip[] are compatible with all routines.
 308 */
 309
 310
 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 }
 329
 330
 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 }
 373
 374
 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 }
 429
 430
 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 }
 485
 486
 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 }
 578
 579
 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 }
 662
 663
 664 /* -------- initializing routines -------- */
 665
 666
 667 #include <math.h>
 668
 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 }
 728
 729
 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 }
 748
 749
 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 }
 767
 768
 769 /* -------- child routines -------- */
 770
 771
 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 <stdio.h>
 782 #define cdft_thread_t pthread_t
 783 #define cdft_thread_create(thp,func,argp) { \
 784     if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
 785         fprintf(stderr, "cdft thread error\n"); \
 786         exit(1); \
 787     } \
 788 }
 789 #define cdft_thread_wait(th) { \
 790     if (pthread_join(th, NULL) != 0) { \
 791         fprintf(stderr, "cdft thread error\n"); \
 792         exit(1); \
 793     } \
 794 }
 795 #endif /* USE_CDFT_PTHREADS */
 796
 797
 798 #ifdef USE_CDFT_WINTHREADS
 799 #define USE_CDFT_THREADS
 800 #ifndef CDFT_THREADS_BEGIN_N
 801 #define CDFT_THREADS_BEGIN_N 32768
 802 #endif
 803 #ifndef CDFT_4THREADS_BEGIN_N
 804 #define CDFT_4THREADS_BEGIN_N 524288
 805 #endif
 806 #include <windows.h>
 807 #include <stdio.h>
 808 #define cdft_thread_t HANDLE
 809 #define cdft_thread_create(thp,func,argp) { \
 810     DWORD thid; \
 811     *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
 812     if (*(thp) == 0) { \
 813         fprintf(stderr, "cdft thread error\n"); \
 814         exit(1); \
 815     } \
 816 }
 817 #define cdft_thread_wait(th) { \
 818     WaitForSingleObject(th, INFINITE); \
 819     CloseHandle(th); \
 820 }
 821 #endif /* USE_CDFT_WINTHREADS */
 822
 823
 824 void cftfsub(int n, double *a, int *ip, int nw, double *w)
 825 {
 826     void bitrv2(int n, int *ip, double *a);
 827     void bitrv216(double *a);
 828     void bitrv208(double *a);
 829     void cftf1st(int n, double *a, double *w);
 830     void cftrec4(int n, double *a, int nw, double *w);
 831     void cftleaf(int n, int isplt, double *a, int nw, double *w);
 832     void cftfx41(int n, double *a, int nw, double *w);
 833     void cftf161(double *a, double *w);
 834     void cftf081(double *a, double *w);
 835     void cftf040(double *a);
 836     void cftx020(double *a);
 837 #ifdef USE_CDFT_THREADS
 838     void cftrec4_th(int n, double *a, int nw, double *w);
 839 #endif /* USE_CDFT_THREADS */
 840
 841     if (n > 8) {
 842         if (n > 32) {
 843             cftf1st(n, a, &w[nw - (n >> 2)]);
 844 #ifdef USE_CDFT_THREADS
 845             if (n > CDFT_THREADS_BEGIN_N) {
 846                 cftrec4_th(n, a, nw, w);
 847             } else
 848 #endif /* USE_CDFT_THREADS */
 849             if (n > 512) {
 850                 cftrec4(n, a, nw, w);
 851             } else if (n > 128) {
 852                 cftleaf(n, 1, a, nw, w);
 853             } else {
 854                 cftfx41(n, a, nw, w);
 855             }
 856             bitrv2(n, ip, a);
 857         } else if (n == 32) {
 858             cftf161(a, &w[nw - 8]);
 859             bitrv216(a);
 860         } else {
 861             cftf081(a, w);
 862             bitrv208(a);
 863         }
 864     } else if (n == 8) {
 865         cftf040(a);
 866     } else if (n == 4) {
 867         cftx020(a);
 868     }
 869 }
 870
 871
 872 void cftbsub(int n, double *a, int *ip, int nw, double *w)
 873 {
 874     void bitrv2conj(int n, int *ip, double *a);
 875     void bitrv216neg(double *a);
 876     void bitrv208neg(double *a);
 877     void cftb1st(int n, double *a, double *w);
 878     void cftrec4(int n, double *a, int nw, double *w);
 879     void cftleaf(int n, int isplt, double *a, int nw, double *w);
 880     void cftfx41(int n, double *a, int nw, double *w);
 881     void cftf161(double *a, double *w);
 882     void cftf081(double *a, double *w);
 883     void cftb040(double *a);
 884     void cftx020(double *a);
 885 #ifdef USE_CDFT_THREADS
 886     void cftrec4_th(int n, double *a, int nw, double *w);
 887 #endif /* USE_CDFT_THREADS */
 888
 889     if (n > 8) {
 890         if (n > 32) {
 891             cftb1st(n, a, &w[nw - (n >> 2)]);
 892 #ifdef USE_CDFT_THREADS
 893             if (n > CDFT_THREADS_BEGIN_N) {
 894                 cftrec4_th(n, a, nw, w);
 895             } else
 896 #endif /* USE_CDFT_THREADS */
 897             if (n > 512) {
 898                 cftrec4(n, a, nw, w);
 899             } else if (n > 128) {
 900                 cftleaf(n, 1, a, nw, w);
 901             } else {
 902                 cftfx41(n, a, nw, w);
 903             }
 904             bitrv2conj(n, ip, a);
 905         } else if (n == 32) {
 906             cftf161(a, &w[nw - 8]);
 907             bitrv216neg(a);
 908         } else {
 909             cftf081(a, w);
 910             bitrv208neg(a);
 911         }
 912     } else if (n == 8) {
 913         cftb040(a);
 914     } else if (n == 4) {
 915         cftx020(a);
 916     }
 917 }
 918
 919
 920 void bitrv2(int n, int *ip, double *a)
 921 {
 922     int j, j1, k, k1, l, m, nh, nm;
 923     double xr, xi, yr, yi;
 924
 925     m = 1;
 926     for (l = n >> 2; l > 8; l >>= 2) {
 927         m <<= 1;
 928     }
 929     nh = n >> 1;
 930     nm = 4 * m;
 931     if (l == 8) {
 932         for (k = 0; k < m; k++) {
 933             for (j = 0; j < k; j++) {
 934                 j1 = 4 * j + 2 * ip[m + k];
 935                 k1 = 4 * k + 2 * ip[m + j];
 936                 xr = a[j1];
 937                 xi = a[j1 + 1];
 938                 yr = a[k1];
 939                 yi = a[k1 + 1];
 940                 a[j1] = yr;
 941                 a[j1 + 1] = yi;
 942                 a[k1] = xr;
 943                 a[k1 + 1] = xi;
 944                 j1 += nm;
 945                 k1 += 2 * nm;
 946                 xr = a[j1];
 947                 xi = a[j1 + 1];
 948                 yr = a[k1];
 949                 yi = a[k1 + 1];
 950                 a[j1] = yr;
 951                 a[j1 + 1] = yi;
 952                 a[k1] = xr;
 953                 a[k1 + 1] = xi;
 954                 j1 += nm;
 955                 k1 -= nm;
 956                 xr = a[j1];
 957                 xi = a[j1 + 1];
 958                 yr = a[k1];
 959                 yi = a[k1 + 1];
 960                 a[j1] = yr;
 961                 a[j1 + 1] = yi;
 962                 a[k1] = xr;
 963                 a[k1 + 1] = xi;
 964                 j1 += nm;
 965                 k1 += 2 * nm;
 966                 xr = a[j1];
 967                 xi = a[j1 + 1];
 968                 yr = a[k1];
 969                 yi = a[k1 + 1];
 970                 a[j1] = yr;
 971                 a[j1 + 1] = yi;
 972                 a[k1] = xr;
 973                 a[k1 + 1] = xi;
 974                 j1 += nh;
 975                 k1 += 2;
 976                 xr = a[j1];
 977                 xi = a[j1 + 1];
 978                 yr = a[k1];
 979                 yi = a[k1 + 1];
 980                 a[j1] = yr;
 981                 a[j1 + 1] = yi;
 982                 a[k1] = xr;
 983                 a[k1 + 1] = xi;
 984                 j1 -= nm;
 985                 k1 -= 2 * nm;
 986                 xr = a[j1];
 987                 xi = a[j1 + 1];
 988                 yr = a[k1];
 989                 yi = a[k1 + 1];
 990                 a[j1] = yr;
 991                 a[j1 + 1] = yi;
 992                 a[k1] = xr;
 993                 a[k1 + 1] = xi;
 994                 j1 -= nm;
 995                 k1 += nm;
 996                 xr = a[j1];
 997                 xi = a[j1 + 1];
 998                 yr = a[k1];
 999                 yi = a[k1 + 1];
1000                 a[j1] = yr;
1001                 a[j1 + 1] = yi;
1002                 a[k1] = xr;
1003                 a[k1 + 1] = xi;
1004                 j1 -= nm;
1005                 k1 -= 2 * nm;
1006                 xr = a[j1];
1007                 xi = a[j1 + 1];
1008                 yr = a[k1];
1009                 yi = a[k1 + 1];
1010                 a[j1] = yr;
1011                 a[j1 + 1] = yi;
1012                 a[k1] = xr;
1013                 a[k1 + 1] = xi;
1014                 j1 += 2;
1015                 k1 += nh;
1016                 xr = a[j1];
1017                 xi = a[j1 + 1];
1018                 yr = a[k1];
1019                 yi = a[k1 + 1];
1020                 a[j1] = yr;
1021                 a[j1 + 1] = yi;
1022                 a[k1] = xr;
1023                 a[k1 + 1] = xi;
1024                 j1 += nm;
1025                 k1 += 2 * nm;
1026                 xr = a[j1];
1027                 xi = a[j1 + 1];
1028                 yr = a[k1];
1029                 yi = a[k1 + 1];
1030                 a[j1] = yr;
1031                 a[j1 + 1] = yi;
1032                 a[k1] = xr;
1033                 a[k1 + 1] = xi;
1034                 j1 += nm;
1035                 k1 -= nm;
1036                 xr = a[j1];
1037                 xi = a[j1 + 1];
1038                 yr = a[k1];
1039                 yi = a[k1 + 1];
1040                 a[j1] = yr;
1041                 a[j1 + 1] = yi;
1042                 a[k1] = xr;
1043                 a[k1 + 1] = xi;
1044                 j1 += nm;
1045                 k1 += 2 * nm;
1046                 xr = a[j1];
1047                 xi = a[j1 + 1];
1048                 yr = a[k1];
1049                 yi = a[k1 + 1];
1050                 a[j1] = yr;
1051                 a[j1 + 1] = yi;
1052                 a[k1] = xr;
1053                 a[k1 + 1] = xi;
1054                 j1 -= nh;
1055                 k1 -= 2;
1056                 xr = a[j1];
1057                 xi = a[j1 + 1];
1058                 yr = a[k1];
1059                 yi = a[k1 + 1];
1060                 a[j1] = yr;
1061                 a[j1 + 1] = yi;
1062                 a[k1] = xr;
1063                 a[k1 + 1] = xi;
1064                 j1 -= nm;
1065                 k1 -= 2 * nm;
1066                 xr = a[j1];
1067                 xi = a[j1 + 1];
1068                 yr = a[k1];
1069                 yi = a[k1 + 1];
1070                 a[j1] = yr;
1071                 a[j1 + 1] = yi;
1072                 a[k1] = xr;
1073                 a[k1 + 1] = xi;
1074                 j1 -= nm;
1075                 k1 += nm;
1076                 xr = a[j1];
1077                 xi = a[j1 + 1];
1078                 yr = a[k1];
1079                 yi = a[k1 + 1];
1080                 a[j1] = yr;
1081                 a[j1 + 1] = yi;
1082                 a[k1] = xr;
1083                 a[k1 + 1] = xi;
1084                 j1 -= nm;
1085                 k1 -= 2 * nm;
1086                 xr = a[j1];
1087                 xi = a[j1 + 1];
1088                 yr = a[k1];
1089                 yi = a[k1 + 1];
1090                 a[j1] = yr;
1091                 a[j1 + 1] = yi;
1092                 a[k1] = xr;
1093                 a[k1 + 1] = xi;
1094             }
1095             k1 = 4 * k + 2 * ip[m + k];
1096             j1 = k1 + 2;
1097             k1 += nh;
1098             xr = a[j1];
1099             xi = a[j1 + 1];
1100             yr = a[k1];
1101             yi = a[k1 + 1];
1102             a[j1] = yr;
1103             a[j1 + 1] = yi;
1104             a[k1] = xr;
1105             a[k1 + 1] = xi;
1106             j1 += nm;
1107             k1 += 2 * nm;
1108             xr = a[j1];
1109             xi = a[j1 + 1];
1110             yr = a[k1];
1111             yi = a[k1 + 1];
1112             a[j1] = yr;
1113             a[j1 + 1] = yi;
1114             a[k1] = xr;
1115             a[k1 + 1] = xi;
1116             j1 += nm;
1117             k1 -= nm;
1118             xr = a[j1];
1119             xi = a[j1 + 1];
1120             yr = a[k1];
1121             yi = a[k1 + 1];
1122             a[j1] = yr;
1123             a[j1 + 1] = yi;
1124             a[k1] = xr;
1125             a[k1 + 1] = xi;
1126             j1 -= 2;
1127             k1 -= nh;
1128             xr = a[j1];
1129             xi = a[j1 + 1];
1130             yr = a[k1];
1131             yi = a[k1 + 1];
1132             a[j1] = yr;
1133             a[j1 + 1] = yi;
1134             a[k1] = xr;
1135             a[k1 + 1] = xi;
1136             j1 += nh + 2;
1137             k1 += nh + 2;
1138             xr = a[j1];
1139             xi = a[j1 + 1];
1140             yr = a[k1];
1141             yi = a[k1 + 1];
1142             a[j1] = yr;
1143             a[j1 + 1] = yi;
1144             a[k1] = xr;
1145             a[k1 + 1] = xi;
1146             j1 -= nh - nm;
1147             k1 += 2 * nm - 2;
1148             xr = a[j1];
1149             xi = a[j1 + 1];
1150             yr = a[k1];
1151             yi = a[k1 + 1];
1152             a[j1] = yr;
1153             a[j1 + 1] = yi;
1154             a[k1] = xr;
1155             a[k1 + 1] = xi;
1156         }
1157     } else {
1158         for (k = 0; k < m; k++) {
1159             for (j = 0; j < k; j++) {
1160                 j1 = 4 * j + ip[m + k];
1161                 k1 = 4 * k + ip[m + j];
1162                 xr = a[j1];
1163                 xi = a[j1 + 1];
1164                 yr = a[k1];
1165                 yi = a[k1 + 1];
1166                 a[j1] = yr;
1167                 a[j1 + 1] = yi;
1168                 a[k1] = xr;
1169                 a[k1 + 1] = xi;
1170                 j1 += nm;
1171                 k1 += nm;
1172                 xr = a[j1];
1173                 xi = a[j1 + 1];
1174                 yr = a[k1];
1175                 yi = a[k1 + 1];
1176                 a[j1] = yr;
1177                 a[j1 + 1] = yi;
1178                 a[k1] = xr;
1179                 a[k1 + 1] = xi;
1180                 j1 += nh;
1181                 k1 += 2;
1182                 xr = a[j1];
1183                 xi = a[j1 + 1];
1184                 yr = a[k1];
1185                 yi = a[k1 + 1];
1186                 a[j1] = yr;
1187                 a[j1 + 1] = yi;
1188                 a[k1] = xr;
1189                 a[k1 + 1] = xi;
1190                 j1 -= nm;
1191                 k1 -= nm;
1192                 xr = a[j1];
1193                 xi = a[j1 + 1];
1194                 yr = a[k1];
1195                 yi = a[k1 + 1];
1196                 a[j1] = yr;
1197                 a[j1 + 1] = yi;
1198                 a[k1] = xr;
1199                 a[k1 + 1] = xi;
1200                 j1 += 2;
1201                 k1 += nh;
1202                 xr = a[j1];
1203                 xi = a[j1 + 1];
1204                 yr = a[k1];
1205                 yi = a[k1 + 1];
1206                 a[j1] = yr;
1207                 a[j1 + 1] = yi;
1208                 a[k1] = xr;
1209                 a[k1 + 1] = xi;
1210                 j1 += nm;
1211                 k1 += nm;
1212                 xr = a[j1];
1213                 xi = a[j1 + 1];
1214                 yr = a[k1];
1215                 yi = a[k1 + 1];
1216                 a[j1] = yr;
1217                 a[j1 + 1] = yi;
1218                 a[k1] = xr;
1219                 a[k1 + 1] = xi;
1220                 j1 -= nh;
1221                 k1 -= 2;
1222                 xr = a[j1];
1223                 xi = a[j1 + 1];
1224                 yr = a[k1];
1225                 yi = a[k1 + 1];
1226                 a[j1] = yr;
1227                 a[j1 + 1] = yi;
1228                 a[k1] = xr;
1229                 a[k1 + 1] = xi;
1230                 j1 -= nm;
1231                 k1 -= nm;
1232                 xr = a[j1];
1233                 xi = a[j1 + 1];
1234                 yr = a[k1];
1235                 yi = a[k1 + 1];
1236                 a[j1] = yr;
1237                 a[j1 + 1] = yi;
1238                 a[k1] = xr;
1239                 a[k1 + 1] = xi;
1240             }
1241             k1 = 4 * k + ip[m + k];
1242             j1 = k1 + 2;
1243             k1 += nh;
1244             xr = a[j1];
1245             xi = a[j1 + 1];
1246             yr = a[k1];
1247             yi = a[k1 + 1];
1248             a[j1] = yr;
1249             a[j1 + 1] = yi;
1250             a[k1] = xr;
1251             a[k1 + 1] = xi;
1252             j1 += nm;
1253             k1 += nm;
1254             xr = a[j1];
1255             xi = a[j1 + 1];
1256             yr = a[k1];
1257             yi = a[k1 + 1];
1258             a[j1] = yr;
1259             a[j1 + 1] = yi;
1260             a[k1] = xr;
1261             a[k1 + 1] = xi;
1262         }
1263     }
1264 }
1265
1266
1267 void bitrv2conj(int n, int *ip, double *a)
1268 {
1269     int j, j1, k, k1, l, m, nh, nm;
1270     double xr, xi, yr, yi;
1271
1272     m = 1;
1273     for (l = n >> 2; l > 8; l >>= 2) {
1274         m <<= 1;
1275     }
1276     nh = n >> 1;
1277     nm = 4 * m;
1278     if (l == 8) {
1279         for (k = 0; k < m; k++) {
1280             for (j = 0; j < k; j++) {
1281                 j1 = 4 * j + 2 * ip[m + k];
1282                 k1 = 4 * k + 2 * ip[m + j];
1283                 xr = a[j1];
1284                 xi = -a[j1 + 1];
1285                 yr = a[k1];
1286                 yi = -a[k1 + 1];
1287                 a[j1] = yr;
1288                 a[j1 + 1] = yi;
1289                 a[k1] = xr;
1290                 a[k1 + 1] = xi;
1291                 j1 += nm;
1292                 k1 += 2 * nm;
1293                 xr = a[j1];
1294                 xi = -a[j1 + 1];
1295                 yr = a[k1];
1296                 yi = -a[k1 + 1];
1297                 a[j1] = yr;
1298                 a[j1 + 1] = yi;
1299                 a[k1] = xr;
1300                 a[k1 + 1] = xi;
1301                 j1 += nm;
1302                 k1 -= nm;
1303                 xr = a[j1];
1304                 xi = -a[j1 + 1];
1305                 yr = a[k1];
1306                 yi = -a[k1 + 1];
1307                 a[j1] = yr;
1308                 a[j1 + 1] = yi;
1309                 a[k1] = xr;
1310                 a[k1 + 1] = xi;
1311                 j1 += nm;
1312                 k1 += 2 * nm;
1313                 xr = a[j1];
1314                 xi = -a[j1 + 1];
1315                 yr = a[k1];
1316                 yi = -a[k1 + 1];
1317                 a[j1] = yr;
1318                 a[j1 + 1] = yi;
1319                 a[k1] = xr;
1320                 a[k1 + 1] = xi;
1321                 j1 += nh;
1322                 k1 += 2;
1323                 xr = a[j1];
1324                 xi = -a[j1 + 1];
1325                 yr = a[k1];
1326                 yi = -a[k1 + 1];
1327                 a[j1] = yr;
1328                 a[j1 + 1] = yi;
1329                 a[k1] = xr;
1330                 a[k1 + 1] = xi;
1331                 j1 -= nm;
1332                 k1 -= 2 * nm;
1333                 xr = a[j1];
1334                 xi = -a[j1 + 1];
1335                 yr = a[k1];
1336                 yi = -a[k1 + 1];
1337                 a[j1] = yr;
1338                 a[j1 + 1] = yi;
1339                 a[k1] = xr;
1340                 a[k1 + 1] = xi;
1341                 j1 -= nm;
1342                 k1 += nm;
1343                 xr = a[j1];
1344                 xi = -a[j1 + 1];
1345                 yr = a[k1];
1346                 yi = -a[k1 + 1];
1347                 a[j1] = yr;
1348                 a[j1 + 1] = yi;
1349                 a[k1] = xr;
1350                 a[k1 + 1] = xi;
1351                 j1 -= nm;
1352                 k1 -= 2 * nm;
1353                 xr = a[j1];
1354                 xi = -a[j1 + 1];
1355                 yr = a[k1];
1356                 yi = -a[k1 + 1];
1357                 a[j1] = yr;
1358                 a[j1 + 1] = yi;
1359                 a[k1] = xr;
1360                 a[k1 + 1] = xi;
1361                 j1 += 2;
1362                 k1 += nh;
1363                 xr = a[j1];
1364                 xi = -a[j1 + 1];
1365                 yr = a[k1];
1366                 yi = -a[k1 + 1];
1367                 a[j1] = yr;
1368                 a[j1 + 1] = yi;
1369                 a[k1] = xr;
1370                 a[k1 + 1] = xi;
1371                 j1 += nm;
1372                 k1 += 2 * nm;
1373                 xr = a[j1];
1374                 xi = -a[j1 + 1];
1375                 yr = a[k1];
1376                 yi = -a[k1 + 1];
1377                 a[j1] = yr;
1378                 a[j1 + 1] = yi;
1379                 a[k1] = xr;
1380                 a[k1 + 1] = xi;
1381                 j1 += nm;
1382                 k1 -= nm;
1383                 xr = a[j1];
1384                 xi = -a[j1 + 1];
1385                 yr = a[k1];
1386                 yi = -a[k1 + 1];
1387                 a[j1] = yr;
1388                 a[j1 + 1] = yi;
1389                 a[k1] = xr;
1390                 a[k1 + 1] = xi;
1391                 j1 += nm;
1392                 k1 += 2 * nm;
1393                 xr = a[j1];
1394                 xi = -a[j1 + 1];
1395                 yr = a[k1];
1396                 yi = -a[k1 + 1];
1397                 a[j1] = yr;
1398                 a[j1 + 1] = yi;
1399                 a[k1] = xr;
1400                 a[k1 + 1] = xi;
1401                 j1 -= nh;
1402                 k1 -= 2;
1403                 xr = a[j1];
1404                 xi = -a[j1 + 1];
1405                 yr = a[k1];
1406                 yi = -a[k1 + 1];
1407                 a[j1] = yr;
1408                 a[j1 + 1] = yi;
1409                 a[k1] = xr;
1410                 a[k1 + 1] = xi;
1411                 j1 -= nm;
1412                 k1 -= 2 * nm;
1413                 xr = a[j1];
1414                 xi = -a[j1 + 1];
1415                 yr = a[k1];
1416                 yi = -a[k1 + 1];
1417                 a[j1] = yr;
1418                 a[j1 + 1] = yi;
1419                 a[k1] = xr;
1420                 a[k1 + 1] = xi;
1421                 j1 -= nm;
1422                 k1 += nm;
1423                 xr = a[j1];
1424                 xi = -a[j1 + 1];
1425                 yr = a[k1];
1426                 yi = -a[k1 + 1];
1427                 a[j1] = yr;
1428                 a[j1 + 1] = yi;
1429                 a[k1] = xr;
1430                 a[k1 + 1] = xi;
1431                 j1 -= nm;
1432                 k1 -= 2 * nm;
1433                 xr = a[j1];
1434                 xi = -a[j1 + 1];
1435                 yr = a[k1];
1436                 yi = -a[k1 + 1];
1437                 a[j1] = yr;
1438                 a[j1 + 1] = yi;
1439                 a[k1] = xr;
1440                 a[k1 + 1] = xi;
1441             }
1442             k1 = 4 * k + 2 * ip[m + k];
1443             j1 = k1 + 2;
1444             k1 += nh;
1445             a[j1 - 1] = -a[j1 - 1];
1446             xr = a[j1];
1447             xi = -a[j1 + 1];
1448             yr = a[k1];
1449             yi = -a[k1 + 1];
1450             a[j1] = yr;
1451             a[j1 + 1] = yi;
1452             a[k1] = xr;
1453             a[k1 + 1] = xi;
1454             a[k1 + 3] = -a[k1 + 3];
1455             j1 += nm;
1456             k1 += 2 * nm;
1457             xr = a[j1];
1458             xi = -a[j1 + 1];
1459             yr = a[k1];
1460             yi = -a[k1 + 1];
1461             a[j1] = yr;
1462             a[j1 + 1] = yi;
1463             a[k1] = xr;
1464             a[k1 + 1] = xi;
1465             j1 += nm;
1466             k1 -= nm;
1467             xr = a[j1];
1468             xi = -a[j1 + 1];
1469             yr = a[k1];
1470             yi = -a[k1 + 1];
1471             a[j1] = yr;
1472             a[j1 + 1] = yi;
1473             a[k1] = xr;
1474             a[k1 + 1] = xi;
1475             j1 -= 2;
1476             k1 -= nh;
1477             xr = a[j1];
1478             xi = -a[j1 + 1];
1479             yr = a[k1];
1480             yi = -a[k1 + 1];
1481             a[j1] = yr;
1482             a[j1 + 1] = yi;
1483             a[k1] = xr;
1484             a[k1 + 1] = xi;
1485             j1 += nh + 2;
1486             k1 += nh + 2;
1487             xr = a[j1];
1488             xi = -a[j1 + 1];
1489             yr = a[k1];
1490             yi = -a[k1 + 1];
1491             a[j1] = yr;
1492             a[j1 + 1] = yi;
1493             a[k1] = xr;
1494             a[k1 + 1] = xi;
1495             j1 -= nh - nm;
1496             k1 += 2 * nm - 2;
1497             a[j1 - 1] = -a[j1 - 1];
1498             xr = a[j1];
1499             xi = -a[j1 + 1];
1500             yr = a[k1];
1501             yi = -a[k1 + 1];
1502             a[j1] = yr;
1503             a[j1 + 1] = yi;
1504             a[k1] = xr;
1505             a[k1 + 1] = xi;
1506             a[k1 + 3] = -a[k1 + 3];
1507         }
1508     } else {
1509         for (k = 0; k < m; k++) {
1510             for (j = 0; j < k; j++) {
1511                 j1 = 4 * j + ip[m + k];
1512                 k1 = 4 * k + ip[m + j];
1513                 xr = a[j1];
1514                 xi = -a[j1 + 1];
1515                 yr = a[k1];
1516                 yi = -a[k1 + 1];
1517                 a[j1] = yr;
1518                 a[j1 + 1] = yi;
1519                 a[k1] = xr;
1520                 a[k1 + 1] = xi;
1521                 j1 += nm;
1522                 k1 += nm;
1523                 xr = a[j1];
1524                 xi = -a[j1 + 1];
1525                 yr = a[k1];
1526                 yi = -a[k1 + 1];
1527                 a[j1] = yr;
1528                 a[j1 + 1] = yi;
1529                 a[k1] = xr;
1530                 a[k1 + 1] = xi;
1531                 j1 += nh;
1532                 k1 += 2;
1533                 xr = a[j1];
1534                 xi = -a[j1 + 1];
1535                 yr = a[k1];
1536                 yi = -a[k1 + 1];
1537                 a[j1] = yr;
1538                 a[j1 + 1] = yi;
1539                 a[k1] = xr;
1540                 a[k1 + 1] = xi;
1541                 j1 -= nm;
1542                 k1 -= nm;
1543                 xr = a[j1];
1544                 xi = -a[j1 + 1];
1545                 yr = a[k1];
1546                 yi = -a[k1 + 1];
1547                 a[j1] = yr;
1548                 a[j1 + 1] = yi;
1549                 a[k1] = xr;
1550                 a[k1 + 1] = xi;
1551                 j1 += 2;
1552                 k1 += nh;
1553                 xr = a[j1];
1554                 xi = -a[j1 + 1];
1555                 yr = a[k1];
1556                 yi = -a[k1 + 1];
1557                 a[j1] = yr;
1558                 a[j1 + 1] = yi;
1559                 a[k1] = xr;
1560                 a[k1 + 1] = xi;
1561                 j1 += nm;
1562                 k1 += nm;
1563                 xr = a[j1];
1564                 xi = -a[j1 + 1];
1565                 yr = a[k1];
1566                 yi = -a[k1 + 1];
1567                 a[j1] = yr;
1568                 a[j1 + 1] = yi;
1569                 a[k1] = xr;
1570                 a[k1 + 1] = xi;
1571                 j1 -= nh;
1572                 k1 -= 2;
1573                 xr = a[j1];
1574                 xi = -a[j1 + 1];
1575                 yr = a[k1];
1576                 yi = -a[k1 + 1];
1577                 a[j1] = yr;
1578                 a[j1 + 1] = yi;
1579                 a[k1] = xr;
1580                 a[k1 + 1] = xi;
1581                 j1 -= nm;
1582                 k1 -= nm;
1583                 xr = a[j1];
1584                 xi = -a[j1 + 1];
1585                 yr = a[k1];
1586                 yi = -a[k1 + 1];
1587                 a[j1] = yr;
1588                 a[j1 + 1] = yi;
1589                 a[k1] = xr;
1590                 a[k1 + 1] = xi;
1591             }
1592             k1 = 4 * k + ip[m + k];
1593             j1 = k1 + 2;
1594             k1 += nh;
1595             a[j1 - 1] = -a[j1 - 1];
1596             xr = a[j1];
1597             xi = -a[j1 + 1];
1598             yr = a[k1];
1599             yi = -a[k1 + 1];
1600             a[j1] = yr;
1601             a[j1 + 1] = yi;
1602             a[k1] = xr;
1603             a[k1 + 1] = xi;
1604             a[k1 + 3] = -a[k1 + 3];
1605             j1 += nm;
1606             k1 += nm;
1607             a[j1 - 1] = -a[j1 - 1];
1608             xr = a[j1];
1609             xi = -a[j1 + 1];
1610             yr = a[k1];
1611             yi = -a[k1 + 1];
1612             a[j1] = yr;
1613             a[j1 + 1] = yi;
1614             a[k1] = xr;
1615             a[k1 + 1] = xi;
1616             a[k1 + 3] = -a[k1 + 3];
1617         }
1618     }
1619 }
1620
1621
1622 void bitrv216(double *a)
1623 {
1624     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1625         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
1626         x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1627
1628     x1r = a[2];
1629     x1i = a[3];
1630     x2r = a[4];
1631     x2i = a[5];
1632     x3r = a[6];
1633     x3i = a[7];
1634     x4r = a[8];
1635     x4i = a[9];
1636     x5r = a[10];
1637     x5i = a[11];
1638     x7r = a[14];
1639     x7i = a[15];
1640     x8r = a[16];
1641     x8i = a[17];
1642     x10r = a[20];
1643     x10i = a[21];
1644     x11r = a[22];
1645     x11i = a[23];
1646     x12r = a[24];
1647     x12i = a[25];
1648     x13r = a[26];
1649     x13i = a[27];
1650     x14r = a[28];
1651     x14i = a[29];
1652     a[2] = x8r;
1653     a[3] = x8i;
1654     a[4] = x4r;
1655     a[5] = x4i;
1656     a[6] = x12r;
1657     a[7] = x12i;
1658     a[8] = x2r;
1659     a[9] = x2i;
1660     a[10] = x10r;
1661     a[11] = x10i;
1662     a[14] = x14r;
1663     a[15] = x14i;
1664     a[16] = x1r;
1665     a[17] = x1i;
1666     a[20] = x5r;
1667     a[21] = x5i;
1668     a[22] = x13r;
1669     a[23] = x13i;
1670     a[24] = x3r;
1671     a[25] = x3i;
1672     a[26] = x11r;
1673     a[27] = x11i;
1674     a[28] = x7r;
1675     a[29] = x7i;
1676 }
1677
1678
1679 void bitrv216neg(double *a)
1680 {
1681     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1682         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
1683         x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
1684         x13r, x13i, x14r, x14i, x15r, x15i;
1685
1686     x1r = a[2];
1687     x1i = a[3];
1688     x2r = a[4];
1689     x2i = a[5];
1690     x3r = a[6];
1691     x3i = a[7];
1692     x4r = a[8];
1693     x4i = a[9];
1694     x5r = a[10];
1695     x5i = a[11];
1696     x6r = a[12];
1697     x6i = a[13];
1698     x7r = a[14];
1699     x7i = a[15];
1700     x8r = a[16];
1701     x8i = a[17];
1702     x9r = a[18];
1703     x9i = a[19];
1704     x10r = a[20];
1705     x10i = a[21];
1706     x11r = a[22];
1707     x11i = a[23];
1708     x12r = a[24];
1709     x12i = a[25];
1710     x13r = a[26];
1711     x13i = a[27];
1712     x14r = a[28];
1713     x14i = a[29];
1714     x15r = a[30];
1715     x15i = a[31];
1716     a[2] = x15r;
1717     a[3] = x15i;
1718     a[4] = x7r;
1719     a[5] = x7i;
1720     a[6] = x11r;
1721     a[7] = x11i;
1722     a[8] = x3r;
1723     a[9] = x3i;
1724     a[10] = x13r;
1725     a[11] = x13i;
1726     a[12] = x5r;
1727     a[13] = x5i;
1728     a[14] = x9r;
1729     a[15] = x9i;
1730     a[16] = x1r;
1731     a[17] = x1i;
1732     a[18] = x14r;
1733     a[19] = x14i;
1734     a[20] = x6r;
1735     a[21] = x6i;
1736     a[22] = x10r;
1737     a[23] = x10i;
1738     a[24] = x2r;
1739     a[25] = x2i;
1740     a[26] = x12r;
1741     a[27] = x12i;
1742     a[28] = x4r;
1743     a[29] = x4i;
1744     a[30] = x8r;
1745     a[31] = x8i;
1746 }
1747
1748
1749 void bitrv208(double *a)
1750 {
1751     double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1752
1753     x1r = a[2];
1754     x1i = a[3];
1755     x3r = a[6];
1756     x3i = a[7];
1757     x4r = a[8];
1758     x4i = a[9];
1759     x6r = a[12];
1760     x6i = a[13];
1761     a[2] = x4r;
1762     a[3] = x4i;
1763     a[6] = x6r;
1764     a[7] = x6i;
1765     a[8] = x1r;
1766     a[9] = x1i;
1767     a[12] = x3r;
1768     a[13] = x3i;
1769 }
1770
1771
1772 void bitrv208neg(double *a)
1773 {
1774     double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1775         x5r, x5i, x6r, x6i, x7r, x7i;
1776
1777     x1r = a[2];
1778     x1i = a[3];
1779     x2r = a[4];
1780     x2i = a[5];
1781     x3r = a[6];
1782     x3i = a[7];
1783     x4r = a[8];
1784     x4i = a[9];
1785     x5r = a[10];
1786     x5i = a[11];
1787     x6r = a[12];
1788     x6i = a[13];
1789     x7r = a[14];
1790     x7i = a[15];
1791     a[2] = x7r;
1792     a[3] = x7i;
1793     a[4] = x3r;
1794     a[5] = x3i;
1795     a[6] = x5r;
1796     a[7] = x5i;
1797     a[8] = x1r;
1798     a[9] = x1i;
1799     a[10] = x6r;
1800     a[11] = x6i;
1801     a[12] = x2r;
1802     a[13] = x2i;
1803     a[14] = x4r;
1804     a[15] = x4i;
1805 }
1806
1807
1808 void cftf1st(int n, double *a, double *w)
1809 {
1810     int j, j0, j1, j2, j3, k, m, mh;
1811     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
1812         wd1r, wd1i, wd3r, wd3i;
1813     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1814         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1815
1816     mh = n >> 3;
1817     m = 2 * mh;
1818     j1 = m;
1819     j2 = j1 + m;
1820     j3 = j2 + m;
1821     x0r = a[0] + a[j2];
1822     x0i = a[1] + a[j2 + 1];
1823     x1r = a[0] - a[j2];
1824     x1i = a[1] - a[j2 + 1];
1825     x2r = a[j1] + a[j3];
1826     x2i = a[j1 + 1] + a[j3 + 1];
1827     x3r = a[j1] - a[j3];
1828     x3i = a[j1 + 1] - a[j3 + 1];
1829     a[0] = x0r + x2r;
1830     a[1] = x0i + x2i;
1831     a[j1] = x0r - x2r;
1832     a[j1 + 1] = x0i - x2i;
1833     a[j2] = x1r - x3i;
1834     a[j2 + 1] = x1i + x3r;
1835     a[j3] = x1r + x3i;
1836     a[j3 + 1] = x1i - x3r;
1837     wn4r = w[1];
1838     csc1 = w[2];
1839     csc3 = w[3];
1840     wd1r = 1;
1841     wd1i = 0;
1842     wd3r = 1;
1843     wd3i = 0;
1844     k = 0;
1845     for (j = 2; j < mh - 2; j += 4) {
1846         k += 4;
1847         wk1r = csc1 * (wd1r + w[k]);
1848         wk1i = csc1 * (wd1i + w[k + 1]);
1849         wk3r = csc3 * (wd3r + w[k + 2]);
1850         wk3i = csc3 * (wd3i + w[k + 3]);
1851         wd1r = w[k];
1852         wd1i = w[k + 1];
1853         wd3r = w[k + 2];
1854         wd3i = w[k + 3];
1855         j1 = j + m;
1856         j2 = j1 + m;
1857         j3 = j2 + m;
1858         x0r = a[j] + a[j2];
1859         x0i = a[j + 1] + a[j2 + 1];
1860         x1r = a[j] - a[j2];
1861         x1i = a[j + 1] - a[j2 + 1];
1862         y0r = a[j + 2] + a[j2 + 2];
1863         y0i = a[j + 3] + a[j2 + 3];
1864         y1r = a[j + 2] - a[j2 + 2];
1865         y1i = a[j + 3] - a[j2 + 3];
1866         x2r = a[j1] + a[j3];
1867         x2i = a[j1 + 1] + a[j3 + 1];
1868         x3r = a[j1] - a[j3];
1869         x3i = a[j1 + 1] - a[j3 + 1];
1870         y2r = a[j1 + 2] + a[j3 + 2];
1871         y2i = a[j1 + 3] + a[j3 + 3];
1872         y3r = a[j1 + 2] - a[j3 + 2];
1873         y3i = a[j1 + 3] - a[j3 + 3];
1874         a[j] = x0r + x2r;
1875         a[j + 1] = x0i + x2i;
1876         a[j + 2] = y0r + y2r;
1877         a[j + 3] = y0i + y2i;
1878         a[j1] = x0r - x2r;
1879         a[j1 + 1] = x0i - x2i;
1880         a[j1 + 2] = y0r - y2r;
1881         a[j1 + 3] = y0i - y2i;
1882         x0r = x1r - x3i;
1883         x0i = x1i + x3r;
1884         a[j2] = wk1r * x0r - wk1i * x0i;
1885         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
1886         x0r = y1r - y3i;
1887         x0i = y1i + y3r;
1888         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
1889         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
1890         x0r = x1r + x3i;
1891         x0i = x1i - x3r;
1892         a[j3] = wk3r * x0r + wk3i * x0i;
1893         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
1894         x0r = y1r + y3i;
1895         x0i = y1i - y3r;
1896         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
1897         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
1898         j0 = m - j;
1899         j1 = j0 + m;
1900         j2 = j1 + m;
1901         j3 = j2 + m;
1902         x0r = a[j0] + a[j2];
1903         x0i = a[j0 + 1] + a[j2 + 1];
1904         x1r = a[j0] - a[j2];
1905         x1i = a[j0 + 1] - a[j2 + 1];
1906         y0r = a[j0 - 2] + a[j2 - 2];
1907         y0i = a[j0 - 1] + a[j2 - 1];
1908         y1r = a[j0 - 2] - a[j2 - 2];
1909         y1i = a[j0 - 1] - a[j2 - 1];
1910         x2r = a[j1] + a[j3];
1911         x2i = a[j1 + 1] + a[j3 + 1];
1912         x3r = a[j1] - a[j3];
1913         x3i = a[j1 + 1] - a[j3 + 1];
1914         y2r = a[j1 - 2] + a[j3 - 2];
1915         y2i = a[j1 - 1] + a[j3 - 1];
1916         y3r = a[j1 - 2] - a[j3 - 2];
1917         y3i = a[j1 - 1] - a[j3 - 1];
1918         a[j0] = x0r + x2r;
1919         a[j0 + 1] = x0i + x2i;
1920         a[j0 - 2] = y0r + y2r;
1921         a[j0 - 1] = y0i + y2i;
1922         a[j1] = x0r - x2r;
1923         a[j1 + 1] = x0i - x2i;
1924         a[j1 - 2] = y0r - y2r;
1925         a[j1 - 1] = y0i - y2i;
1926         x0r = x1r - x3i;
1927         x0i = x1i + x3r;
1928         a[j2] = wk1i * x0r - wk1r * x0i;
1929         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
1930         x0r = y1r - y3i;
1931         x0i = y1i + y3r;
1932         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
1933         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
1934         x0r = x1r + x3i;
1935         x0i = x1i - x3r;
1936         a[j3] = wk3i * x0r + wk3r * x0i;
1937         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
1938         x0r = y1r + y3i;
1939         x0i = y1i - y3r;
1940         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
1941         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
1942     }
1943     wk1r = csc1 * (wd1r + wn4r);
1944     wk1i = csc1 * (wd1i + wn4r);
1945     wk3r = csc3 * (wd3r - wn4r);
1946     wk3i = csc3 * (wd3i - wn4r);
1947     j0 = mh;
1948     j1 = j0 + m;
1949     j2 = j1 + m;
1950     j3 = j2 + m;
1951     x0r = a[j0 - 2] + a[j2 - 2];
1952     x0i = a[j0 - 1] + a[j2 - 1];
1953     x1r = a[j0 - 2] - a[j2 - 2];
1954     x1i = a[j0 - 1] - a[j2 - 1];
1955     x2r = a[j1 - 2] + a[j3 - 2];
1956     x2i = a[j1 - 1] + a[j3 - 1];
1957     x3r = a[j1 - 2] - a[j3 - 2];
1958     x3i = a[j1 - 1] - a[j3 - 1];
1959     a[j0 - 2] = x0r + x2r;
1960     a[j0 - 1] = x0i + x2i;
1961     a[j1 - 2] = x0r - x2r;
1962     a[j1 - 1] = x0i - x2i;
1963     x0r = x1r - x3i;
1964     x0i = x1i + x3r;
1965     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
1966     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
1967     x0r = x1r + x3i;
1968     x0i = x1i - x3r;
1969     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
1970     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
1971     x0r = a[j0] + a[j2];
1972     x0i = a[j0 + 1] + a[j2 + 1];
1973     x1r = a[j0] - a[j2];
1974     x1i = a[j0 + 1] - a[j2 + 1];
1975     x2r = a[j1] + a[j3];
1976     x2i = a[j1 + 1] + a[j3 + 1];
1977     x3r = a[j1] - a[j3];
1978     x3i = a[j1 + 1] - a[j3 + 1];
1979     a[j0] = x0r + x2r;
1980     a[j0 + 1] = x0i + x2i;
1981     a[j1] = x0r - x2r;
1982     a[j1 + 1] = x0i - x2i;
1983     x0r = x1r - x3i;
1984     x0i = x1i + x3r;
1985     a[j2] = wn4r * (x0r - x0i);
1986     a[j2 + 1] = wn4r * (x0i + x0r);
1987     x0r = x1r + x3i;
1988     x0i = x1i - x3r;
1989     a[j3] = -wn4r * (x0r + x0i);
1990     a[j3 + 1] = -wn4r * (x0i - x0r);
1991     x0r = a[j0 + 2] + a[j2 + 2];
1992     x0i = a[j0 + 3] + a[j2 + 3];
1993     x1r = a[j0 + 2] - a[j2 + 2];
1994     x1i = a[j0 + 3] - a[j2 + 3];
1995     x2r = a[j1 + 2] + a[j3 + 2];
1996     x2i = a[j1 + 3] + a[j3 + 3];
1997     x3r = a[j1 + 2] - a[j3 + 2];
1998     x3i = a[j1 + 3] - a[j3 + 3];
1999     a[j0 + 2] = x0r + x2r;
2000     a[j0 + 3] = x0i + x2i;
2001     a[j1 + 2] = x0r - x2r;
2002     a[j1 + 3] = x0i - x2i;
2003     x0r = x1r - x3i;
2004     x0i = x1i + x3r;
2005     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2006     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2007     x0r = x1r + x3i;
2008     x0i = x1i - x3r;
2009     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2010     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2011 }
2012
2013
2014 void cftb1st(int n, double *a, double *w)
2015 {
2016     int j, j0, j1, j2, j3, k, m, mh;
2017     double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
2018         wd1r, wd1i, wd3r, wd3i;
2019     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2020         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
2021
2022     mh = n >> 3;
2023     m = 2 * mh;
2024     j1 = m;
2025     j2 = j1 + m;
2026     j3 = j2 + m;
2027     x0r = a[0] + a[j2];
2028     x0i = -a[1] - a[j2 + 1];
2029     x1r = a[0] - a[j2];
2030     x1i = -a[1] + a[j2 + 1];
2031     x2r = a[j1] + a[j3];
2032     x2i = a[j1 + 1] + a[j3 + 1];
2033     x3r = a[j1] - a[j3];
2034     x3i = a[j1 + 1] - a[j3 + 1];
2035     a[0] = x0r + x2r;
2036     a[1] = x0i - x2i;
2037     a[j1] = x0r - x2r;
2038     a[j1 + 1] = x0i + x2i;
2039     a[j2] = x1r + x3i;
2040     a[j2 + 1] = x1i + x3r;
2041     a[j3] = x1r - x3i;
2042     a[j3 + 1] = x1i - x3r;
2043     wn4r = w[1];
2044     csc1 = w[2];
2045     csc3 = w[3];
2046     wd1r = 1;
2047     wd1i = 0;
2048     wd3r = 1;
2049     wd3i = 0;
2050     k = 0;
2051     for (j = 2; j < mh - 2; j += 4) {
2052         k += 4;
2053         wk1r = csc1 * (wd1r + w[k]);
2054         wk1i = csc1 * (wd1i + w[k + 1]);
2055         wk3r = csc3 * (wd3r + w[k + 2]);
2056         wk3i = csc3 * (wd3i + w[k + 3]);
2057         wd1r = w[k];
2058         wd1i = w[k + 1];
2059         wd3r = w[k + 2];
2060         wd3i = w[k + 3];
2061         j1 = j + m;
2062         j2 = j1 + m;
2063         j3 = j2 + m;
2064         x0r = a[j] + a[j2];
2065         x0i = -a[j + 1] - a[j2 + 1];
2066         x1r = a[j] - a[j2];
2067         x1i = -a[j + 1] + a[j2 + 1];
2068         y0r = a[j + 2] + a[j2 + 2];
2069         y0i = -a[j + 3] - a[j2 + 3];
2070         y1r = a[j + 2] - a[j2 + 2];
2071         y1i = -a[j + 3] + a[j2 + 3];
2072         x2r = a[j1] + a[j3];
2073         x2i = a[j1 + 1] + a[j3 + 1];
2074         x3r = a[j1] - a[j3];
2075         x3i = a[j1 + 1] - a[j3 + 1];
2076         y2r = a[j1 + 2] + a[j3 + 2];
2077         y2i = a[j1 + 3] + a[j3 + 3];
2078         y3r = a[j1 + 2] - a[j3 + 2];
2079         y3i = a[j1 + 3] - a[j3 + 3];
2080         a[j] = x0r + x2r;
2081         a[j + 1] = x0i - x2i;
2082         a[j + 2] = y0r + y2r;
2083         a[j + 3] = y0i - y2i;
2084         a[j1] = x0r - x2r;
2085         a[j1 + 1] = x0i + x2i;
2086         a[j1 + 2] = y0r - y2r;
2087         a[j1 + 3] = y0i + y2i;
2088         x0r = x1r + x3i;
2089         x0i = x1i + x3r;
2090         a[j2] = wk1r * x0r - wk1i * x0i;
2091         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2092         x0r = y1r + y3i;
2093         x0i = y1i + y3r;
2094         a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2095         a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2096         x0r = x1r - x3i;
2097         x0i = x1i - x3r;
2098         a[j3] = wk3r * x0r + wk3i * x0i;
2099         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2100         x0r = y1r - y3i;
2101         x0i = y1i - y3r;
2102         a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2103         a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2104         j0 = m - j;
2105         j1 = j0 + m;
2106         j2 = j1 + m;
2107         j3 = j2 + m;
2108         x0r = a[j0] + a[j2];
2109         x0i = -a[j0 + 1] - a[j2 + 1];
2110         x1r = a[j0] - a[j2];
2111         x1i = -a[j0 + 1] + a[j2 + 1];
2112         y0r = a[j0 - 2] + a[j2 - 2];
2113         y0i = -a[j0 - 1] - a[j2 - 1];
2114         y1r = a[j0 - 2] - a[j2 - 2];
2115         y1i = -a[j0 - 1] + a[j2 - 1];
2116         x2r = a[j1] + a[j3];
2117         x2i = a[j1 + 1] + a[j3 + 1];
2118         x3r = a[j1] - a[j3];
2119         x3i = a[j1 + 1] - a[j3 + 1];
2120         y2r = a[j1 - 2] + a[j3 - 2];
2121         y2i = a[j1 - 1] + a[j3 - 1];
2122         y3r = a[j1 - 2] - a[j3 - 2];
2123         y3i = a[j1 - 1] - a[j3 - 1];
2124         a[j0] = x0r + x2r;
2125         a[j0 + 1] = x0i - x2i;
2126         a[j0 - 2] = y0r + y2r;
2127         a[j0 - 1] = y0i - y2i;
2128         a[j1] = x0r - x2r;
2129         a[j1 + 1] = x0i + x2i;
2130         a[j1 - 2] = y0r - y2r;
2131         a[j1 - 1] = y0i + y2i;
2132         x0r = x1r + x3i;
2133         x0i = x1i + x3r;
2134         a[j2] = wk1i * x0r - wk1r * x0i;
2135         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2136         x0r = y1r + y3i;
2137         x0i = y1i + y3r;
2138         a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2139         a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2140         x0r = x1r - x3i;
2141         x0i = x1i - x3r;
2142         a[j3] = wk3i * x0r + wk3r * x0i;
2143         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2144         x0r = y1r - y3i;
2145         x0i = y1i - y3r;
2146         a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2147         a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2148     }
2149     wk1r = csc1 * (wd1r + wn4r);
2150     wk1i = csc1 * (wd1i + wn4r);
2151     wk3r = csc3 * (wd3r - wn4r);
2152     wk3i = csc3 * (wd3i - wn4r);
2153     j0 = mh;
2154     j1 = j0 + m;
2155     j2 = j1 + m;
2156     j3 = j2 + m;
2157     x0r = a[j0 - 2] + a[j2 - 2];
2158     x0i = -a[j0 - 1] - a[j2 - 1];
2159     x1r = a[j0 - 2] - a[j2 - 2];
2160     x1i = -a[j0 - 1] + a[j2 - 1];
2161     x2r = a[j1 - 2] + a[j3 - 2];
2162     x2i = a[j1 - 1] + a[j3 - 1];
2163     x3r = a[j1 - 2] - a[j3 - 2];
2164     x3i = a[j1 - 1] - a[j3 - 1];
2165     a[j0 - 2] = x0r + x2r;
2166     a[j0 - 1] = x0i - x2i;
2167     a[j1 - 2] = x0r - x2r;
2168     a[j1 - 1] = x0i + x2i;
2169     x0r = x1r + x3i;
2170     x0i = x1i + x3r;
2171     a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2172     a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2173     x0r = x1r - x3i;
2174     x0i = x1i - x3r;
2175     a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2176     a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2177     x0r = a[j0] + a[j2];
2178     x0i = -a[j0 + 1] - a[j2 + 1];
2179     x1r = a[j0] - a[j2];
2180     x1i = -a[j0 + 1] + a[j2 + 1];
2181     x2r = a[j1] + a[j3];
2182     x2i = a[j1 + 1] + a[j3 + 1];
2183     x3r = a[j1] - a[j3];
2184     x3i = a[j1 + 1] - a[j3 + 1];
2185     a[j0] = x0r + x2r;
2186     a[j0 + 1] = x0i - x2i;
2187     a[j1] = x0r - x2r;
2188     a[j1 + 1] = x0i + x2i;
2189     x0r = x1r + x3i;
2190     x0i = x1i + x3r;
2191     a[j2] = wn4r * (x0r - x0i);
2192     a[j2 + 1] = wn4r * (x0i + x0r);
2193     x0r = x1r - x3i;
2194     x0i = x1i - x3r;
2195     a[j3] = -wn4r * (x0r + x0i);
2196     a[j3 + 1] = -wn4r * (x0i - x0r);
2197     x0r = a[j0 + 2] + a[j2 + 2];
2198     x0i = -a[j0 + 3] - a[j2 + 3];
2199     x1r = a[j0 + 2] - a[j2 + 2];
2200     x1i = -a[j0 + 3] + a[j2 + 3];
2201     x2r = a[j1 + 2] + a[j3 + 2];
2202     x2i = a[j1 + 3] + a[j3 + 3];
2203     x3r = a[j1 + 2] - a[j3 + 2];
2204     x3i = a[j1 + 3] - a[j3 + 3];
2205     a[j0 + 2] = x0r + x2r;
2206     a[j0 + 3] = x0i - x2i;
2207     a[j1 + 2] = x0r - x2r;
2208     a[j1 + 3] = x0i + x2i;
2209     x0r = x1r + x3i;
2210     x0i = x1i + x3r;
2211     a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2212     a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2213     x0r = x1r - x3i;
2214     x0i = x1i - x3r;
2215     a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2216     a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2217 }
2218
2219
2220 #ifdef USE_CDFT_THREADS
2221 struct cdft_arg_st {
2222     int n0;
2223     int n;
2224     double *a;
2225     int nw;
2226     double *w;
2227 };
2228 typedef struct cdft_arg_st cdft_arg_t;
2229
2230
2231 void cftrec4_th(int n, double *a, int nw, double *w)
2232 {
2233     void *cftrec1_th(void *p);
2234     void *cftrec2_th(void *p);
2235     int i, idiv4, m, nthread;
2236     cdft_thread_t th[4];
2237     cdft_arg_t ag[4];
2238
2239     nthread = 2;
2240     idiv4 = 0;
2241     m = n >> 1;
2242     if (n > CDFT_4THREADS_BEGIN_N) {
2243         nthread = 4;
2244         idiv4 = 1;
2245         m >>= 1;
2246     }
2247     for (i = 0; i < nthread; i++) {
2248         ag[i].n0 = n;
2249         ag[i].n = m;
2250         ag[i].a = &a[i * m];
2251         ag[i].nw = nw;
2252         ag[i].w = w;
2253         if (i != idiv4) {
2254             cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
2255         } else {
2256             cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
2257         }
2258     }
2259     for (i = 0; i < nthread; i++) {
2260         cdft_thread_wait(th[i]);
2261     }
2262 }
2263
2264
2265 void *cftrec1_th(void *p)
2266 {
2267     int cfttree(int n, int j, int k, double *a, int nw, double *w);
2268     void cftleaf(int n, int isplt, double *a, int nw, double *w);
2269     void cftmdl1(int n, double *a, double *w);
2270     int isplt, j, k, m, n, n0, nw;
2271     double *a, *w;
2272
2273     n0 = ((cdft_arg_t *) p)->n0;
2274     n = ((cdft_arg_t *) p)->n;
2275     a = ((cdft_arg_t *) p)->a;
2276     nw = ((cdft_arg_t *) p)->nw;
2277     w = ((cdft_arg_t *) p)->w;
2278     m = n0;
2279     while (m > 512) {
2280         m >>= 2;
2281         cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2282     }
2283     cftleaf(m, 1, &a[n - m], nw, w);
2284     k = 0;
2285     for (j = n - m; j > 0; j -= m) {
2286         k++;
2287         isplt = cfttree(m, j, k, a, nw, w);
2288         cftleaf(m, isplt, &a[j - m], nw, w);
2289     }
2290     return (void *) 0;
2291 }
2292
2293
2294 void *cftrec2_th(void *p)
2295 {
2296     int cfttree(int n, int j, int k, double *a, int nw, double *w);
2297     void cftleaf(int n, int isplt, double *a, int nw, double *w);
2298     void cftmdl2(int n, double *a, double *w);
2299     int isplt, j, k, m, n, n0, nw;
2300     double *a, *w;
2301
2302     n0 = ((cdft_arg_t *) p)->n0;
2303     n = ((cdft_arg_t *) p)->n;
2304     a = ((cdft_arg_t *) p)->a;
2305     nw = ((cdft_arg_t *) p)->nw;
2306     w = ((cdft_arg_t *) p)->w;
2307     k = 1;
2308     m = n0;
2309     while (m > 512) {
2310         m >>= 2;
2311         k <<= 2;
2312         cftmdl2(m, &a[n - m], &w[nw - m]);
2313     }
2314     cftleaf(m, 0, &a[n - m], nw, w);
2315     k >>= 1;
2316     for (j = n - m; j > 0; j -= m) {
2317         k++;
2318         isplt = cfttree(m, j, k, a, nw, w);
2319         cftleaf(m, isplt, &a[j - m], nw, w);
2320     }
2321     return (void *) 0;
2322 }
2323 #endif /* USE_CDFT_THREADS */
2324
2325
2326 void cftrec4(int n, double *a, int nw, double *w)
2327 {
2328     int cfttree(int n, int j, int k, double *a, int nw, double *w);
2329     void cftleaf(int n, int isplt, double *a, int nw, double *w);
2330     void cftmdl1(int n, double *a, double *w);
2331     int isplt, j, k, m;
2332
2333     m = n;
2334     while (m > 512) {
2335         m >>= 2;
2336         cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2337     }
2338     cftleaf(m, 1, &a[n - m], nw, w);
2339     k = 0;
2340     for (j = n - m; j > 0; j -= m) {
2341         k++;
2342         isplt = cfttree(m, j, k, a, nw, w);
2343         cftleaf(m, isplt, &a[j - m], nw, w);
2344     }
2345 }
2346
2347
2348 int cfttree(int n, int j, int k, double *a, int nw, double *w)
2349 {
2350     void cftmdl1(int n, double *a, double *w);
2351     void cftmdl2(int n, double *a, double *w);
2352     int i, isplt, m;
2353
2354     if ((k & 3) != 0) {
2355         isplt = k & 1;
2356         if (isplt != 0) {
2357             cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
2358         } else {
2359             cftmdl2(n, &a[j - n], &w[nw - n]);
2360         }
2361     } else {
2362         m = n;
2363         for (i = k; (i & 3) == 0; i >>= 2) {
2364             m <<= 2;
2365         }
2366         isplt = i & 1;
2367         if (isplt != 0) {
2368             while (m > 128) {
2369                 cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
2370                 m >>= 2;
2371             }
2372         } else {
2373             while (m > 128) {
2374                 cftmdl2(m, &a[j - m], &w[nw - m]);
2375                 m >>= 2;
2376             }
2377         }
2378     }
2379     return isplt;
2380 }
2381
2382
2383 void cftleaf(int n, int isplt, double *a, int nw, double *w)
2384 {
2385     void cftmdl1(int n, double *a, double *w);
2386     void cftmdl2(int n, double *a, double *w);
2387     void cftf161(double *a, double *w);
2388     void cftf162(double *a, double *w);
2389     void cftf081(double *a, double *w);
2390     void cftf082(double *a, double *w);
2391
2392     if (n == 512) {
2393         cftmdl1(128, a, &w[nw - 64]);
2394         cftf161(a, &w[nw - 8]);
2395         cftf162(&a[32], &w[nw - 32]);
2396         cftf161(&a[64], &w[nw - 8]);
2397         cftf161(&a[96], &w[nw - 8]);
2398         cftmdl2(128, &a[128], &w[nw - 128]);
2399         cftf161(&a[128], &w[nw - 8]);
2400         cftf162(&a[160], &w[nw - 32]);
2401         cftf161(&a[192], &w[nw - 8]);
2402         cftf162(&a[224], &w[nw - 32]);
2403         cftmdl1(128, &a[256], &w[nw - 64]);
2404         cftf161(&a[256], &w[nw - 8]);
2405         cftf162(&a[288], &w[nw - 32]);
2406         cftf161(&a[320], &w[nw - 8]);
2407         cftf161(&a[352], &w[nw - 8]);
2408         if (isplt != 0) {
2409             cftmdl1(128, &a[384], &w[nw - 64]);
2410             cftf161(&a[480], &w[nw - 8]);
2411         } else {
2412             cftmdl2(128, &a[384], &w[nw - 128]);
2413             cftf162(&a[480], &w[nw - 32]);
2414         }
2415         cftf161(&a[384], &w[nw - 8]);
2416         cftf162(&a[416], &w[nw - 32]);
2417         cftf161(&a[448], &w[nw - 8]);
2418     } else {
2419         cftmdl1(64, a, &w[nw - 32]);
2420         cftf081(a, &w[nw - 8]);
2421         cftf082(&a[16], &w[nw - 8]);
2422         cftf081(&a[32], &w[nw - 8]);
2423         cftf081(&a[48], &w[nw - 8]);
2424         cftmdl2(64, &a[64], &w[nw - 64]);
2425         cftf081(&a[64], &w[nw - 8]);
2426         cftf082(&a[80], &w[nw - 8]);
2427         cftf081(&a[96], &w[nw - 8]);
2428         cftf082(&a[112], &w[nw - 8]);
2429         cftmdl1(64, &a[128], &w[nw - 32]);
2430         cftf081(&a[128], &w[nw - 8]);
2431         cftf082(&a[144], &w[nw - 8]);
2432         cftf081(&a[160], &w[nw - 8]);
2433         cftf081(&a[176], &w[nw - 8]);
2434         if (isplt != 0) {
2435             cftmdl1(64, &a[192], &w[nw - 32]);
2436             cftf081(&a[240], &w[nw - 8]);
2437         } else {
2438             cftmdl2(64, &a[192], &w[nw - 64]);
2439             cftf082(&a[240], &w[nw - 8]);
2440         }
2441         cftf081(&a[192], &w[nw - 8]);
2442         cftf082(&a[208], &w[nw - 8]);
2443         cftf081(&a[224], &w[nw - 8]);
2444     }
2445 }
2446
2447
2448 void cftmdl1(int n, double *a, double *w)
2449 {
2450     int j, j0, j1, j2, j3, k, m, mh;
2451     double wn4r, wk1r, wk1i, wk3r, wk3i;
2452     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2453
2454     mh = n >> 3;
2455     m = 2 * mh;
2456     j1 = m;
2457     j2 = j1 + m;
2458     j3 = j2 + m;
2459     x0r = a[0] + a[j2];
2460     x0i = a[1] + a[j2 + 1];
2461     x1r = a[0] - a[j2];
2462     x1i = a[1] - a[j2 + 1];
2463     x2r = a[j1] + a[j3];
2464     x2i = a[j1 + 1] + a[j3 + 1];
2465     x3r = a[j1] - a[j3];
2466     x3i = a[j1 + 1] - a[j3 + 1];
2467     a[0] = x0r + x2r;
2468     a[1] = x0i + x2i;
2469     a[j1] = x0r - x2r;
2470     a[j1 + 1] = x0i - x2i;
2471     a[j2] = x1r - x3i;
2472     a[j2 + 1] = x1i + x3r;
2473     a[j3] = x1r + x3i;
2474     a[j3 + 1] = x1i - x3r;
2475     wn4r = w[1];
2476     k = 0;
2477     for (j = 2; j < mh; j += 2) {
2478         k += 4;
2479         wk1r = w[k];
2480         wk1i = w[k + 1];
2481         wk3r = w[k + 2];
2482         wk3i = w[k + 3];
2483         j1 = j + m;
2484         j2 = j1 + m;
2485         j3 = j2 + m;
2486         x0r = a[j] + a[j2];
2487         x0i = a[j + 1] + a[j2 + 1];
2488         x1r = a[j] - a[j2];
2489         x1i = a[j + 1] - a[j2 + 1];
2490         x2r = a[j1] + a[j3];
2491         x2i = a[j1 + 1] + a[j3 + 1];
2492         x3r = a[j1] - a[j3];
2493         x3i = a[j1 + 1] - a[j3 + 1];
2494         a[j] = x0r + x2r;
2495         a[j + 1] = x0i + x2i;
2496         a[j1] = x0r - x2r;
2497         a[j1 + 1] = x0i - x2i;
2498         x0r = x1r - x3i;
2499         x0i = x1i + x3r;
2500         a[j2] = wk1r * x0r - wk1i * x0i;
2501         a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2502         x0r = x1r + x3i;
2503         x0i = x1i - x3r;
2504         a[j3] = wk3r * x0r + wk3i * x0i;
2505         a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2506         j0 = m - j;
2507         j1 = j0 + m;
2508         j2 = j1 + m;
2509         j3 = j2 + m;
2510         x0r = a[j0] + a[j2];
2511         x0i = a[j0 + 1] + a[j2 + 1];
2512         x1r = a[j0] - a[j2];
2513         x1i = a[j0 + 1] - a[j2 + 1];
2514         x2r = a[j1] + a[j3];
2515         x2i = a[j1 + 1] + a[j3 + 1];
2516         x3r = a[j1] - a[j3];
2517         x3i = a[j1 + 1] - a[j3 + 1];
2518         a[j0] = x0r + x2r;
2519         a[j0 + 1] = x0i + x2i;
2520         a[j1] = x0r - x2r;
2521         a[j1 + 1] = x0i - x2i;
2522         x0r = x1r - x3i;
2523         x0i = x1i + x3r;
2524         a[j2] = wk1i * x0r - wk1r * x0i;
2525         a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2526         x0r = x1r + x3i;
2527         x0i = x1i - x3r;
2528         a[j3] = wk3i * x0r + wk3r * x0i;
2529         a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2530     }
2531     j0 = mh;
2532     j1 = j0 + m;
2533     j2 = j1 + m;
2534     j3 = j2 + m;
2535     x0r = a[j0] + a[j2];
2536     x0i = a[j0 + 1] + a[j2 + 1];
2537     x1r = a[j0] - a[j2];
2538     x1i = a[j0 + 1] - a[j2 + 1];
2539     x2r = a[j1] + a[j3];
2540     x2i = a[j1 + 1] + a[j3 + 1];
2541     x3r = a[j1] - a[j3];
2542     x3i = a[j1 + 1] - a[j3 + 1];
2543     a[j0] = x0r + x2r;
2544     a[j0 + 1] = x0i + x2i;
2545     a[j1] = x0r - x2r;
2546     a[j1 + 1] = x0i - x2i;
2547     x0r = x1r - x3i;
2548     x0i = x1i + x3r;
2549     a[j2] = wn4r * (x0r - x0i);
2550     a[j2 + 1] = wn4r * (x0i + x0r);
2551     x0r = x1r + x3i;
2552     x0i = x1i - x3r;
2553     a[j3] = -wn4r * (x0r + x0i);
2554     a[j3 + 1] = -wn4r * (x0i - x0r);
2555 }
2556
2557
2558 void cftmdl2(int n, double *a, double *w)
2559 {
2560     int j, j0, j1, j2, j3, k, kr, m, mh;
2561     double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2562     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2563
2564     mh = n >> 3;
2565     m = 2 * mh;
2566     wn4r = w[1];
2567     j1 = m;
2568     j2 = j1 + m;
2569     j3 = j2 + m;
2570     x0r = a[0] - a[j2 + 1];
2571     x0i = a[1] + a[j2];
2572     x1r = a[0] + a[j2 + 1];
2573     x1i = a[1] - a[j2];
2574     x2r = a[j1] - a[j3 + 1];
2575     x2i = a[j1 + 1] + a[j3];
2576     x3r = a[j1] + a[j3 + 1];
2577     x3i = a[j1 + 1] - a[j3];
2578     y0r = wn4r * (x2r - x2i);
2579     y0i = wn4r * (x2i + x2r);
2580     a[0] = x0r + y0r;
2581     a[1] = x0i + y0i;
2582     a[j1] = x0r - y0r;
2583     a[j1 + 1] = x0i - y0i;
2584     y0r = wn4r * (x3r - x3i);
2585     y0i = wn4r * (x3i + x3r);
2586     a[j2] = x1r - y0i;
2587     a[j2 + 1] = x1i + y0r;
2588     a[j3] = x1r + y0i;
2589     a[j3 + 1] = x1i - y0r;
2590     k = 0;
2591     kr = 2 * m;
2592     for (j = 2; j < mh; j += 2) {
2593         k += 4;
2594         wk1r = w[k];
2595         wk1i = w[k + 1];
2596         wk3r = w[k + 2];
2597         wk3i = w[k + 3];
2598         kr -= 4;
2599         wd1i = w[kr];
2600         wd1r = w[kr + 1];
2601         wd3i = w[kr + 2];
2602         wd3r = w[kr + 3];
2603         j1 = j + m;
2604         j2 = j1 + m;
2605         j3 = j2 + m;
2606         x0r = a[j] - a[j2 + 1];
2607         x0i = a[j + 1] + a[j2];
2608         x1r = a[j] + a[j2 + 1];
2609         x1i = a[j + 1] - a[j2];
2610         x2r = a[j1] - a[j3 + 1];
2611         x2i = a[j1 + 1] + a[j3];
2612         x3r = a[j1] + a[j3 + 1];
2613         x3i = a[j1 + 1] - a[j3];
2614         y0r = wk1r * x0r - wk1i * x0i;
2615         y0i = wk1r * x0i + wk1i * x0r;
2616         y2r = wd1r * x2r - wd1i * x2i;
2617         y2i = wd1r * x2i + wd1i * x2r;
2618         a[j] = y0r + y2r;
2619         a[j + 1] = y0i + y2i;
2620         a[j1] = y0r - y2r;
2621         a[j1 + 1] = y0i - y2i;
2622         y0r = wk3r * x1r + wk3i * x1i;
2623         y0i = wk3r * x1i - wk3i * x1r;
2624         y2r = wd3r * x3r + wd3i * x3i;
2625         y2i = wd3r * x3i - wd3i * x3r;
2626         a[j2] = y0r + y2r;
2627         a[j2 + 1] = y0i + y2i;
2628         a[j3] = y0r - y2r;
2629         a[j3 + 1] = y0i - y2i;
2630         j0 = m - j;
2631         j1 = j0 + m;
2632         j2 = j1 + m;
2633         j3 = j2 + m;
2634         x0r = a[j0] - a[j2 + 1];
2635         x0i = a[j0 + 1] + a[j2];
2636         x1r = a[j0] + a[j2 + 1];
2637         x1i = a[j0 + 1] - a[j2];
2638         x2r = a[j1] - a[j3 + 1];
2639         x2i = a[j1 + 1] + a[j3];
2640         x3r = a[j1] + a[j3 + 1];
2641         x3i = a[j1 + 1] - a[j3];
2642         y0r = wd1i * x0r - wd1r * x0i;
2643         y0i = wd1i * x0i + wd1r * x0r;
2644         y2r = wk1i * x2r - wk1r * x2i;
2645         y2i = wk1i * x2i + wk1r * x2r;
2646         a[j0] = y0r + y2r;
2647         a[j0 + 1] = y0i + y2i;
2648         a[j1] = y0r - y2r;
2649         a[j1 + 1] = y0i - y2i;
2650         y0r = wd3i * x1r + wd3r * x1i;
2651         y0i = wd3i * x1i - wd3r * x1r;
2652         y2r = wk3i * x3r + wk3r * x3i;
2653         y2i = wk3i * x3i - wk3r * x3r;
2654         a[j2] = y0r + y2r;
2655         a[j2 + 1] = y0i + y2i;
2656         a[j3] = y0r - y2r;
2657         a[j3 + 1] = y0i - y2i;
2658     }
2659     wk1r = w[m];
2660     wk1i = w[m + 1];
2661     j0 = mh;
2662     j1 = j0 + m;
2663     j2 = j1 + m;
2664     j3 = j2 + m;
2665     x0r = a[j0] - a[j2 + 1];
2666     x0i = a[j0 + 1] + a[j2];
2667     x1r = a[j0] + a[j2 + 1];
2668     x1i = a[j0 + 1] - a[j2];
2669     x2r = a[j1] - a[j3 + 1];
2670     x2i = a[j1 + 1] + a[j3];
2671     x3r = a[j1] + a[j3 + 1];
2672     x3i = a[j1 + 1] - a[j3];
2673     y0r = wk1r * x0r - wk1i * x0i;
2674     y0i = wk1r * x0i + wk1i * x0r;
2675     y2r = wk1i * x2r - wk1r * x2i;
2676     y2i = wk1i * x2i + wk1r * x2r;
2677     a[j0] = y0r + y2r;
2678     a[j0 + 1] = y0i + y2i;
2679     a[j1] = y0r - y2r;
2680     a[j1 + 1] = y0i - y2i;
2681     y0r = wk1i * x1r - wk1r * x1i;
2682     y0i = wk1i * x1i + wk1r * x1r;
2683     y2r = wk1r * x3r - wk1i * x3i;
2684     y2i = wk1r * x3i + wk1i * x3r;
2685     a[j2] = y0r - y2r;
2686     a[j2 + 1] = y0i - y2i;
2687     a[j3] = y0r + y2r;
2688     a[j3 + 1] = y0i + y2i;
2689 }
2690
2691
2692 void cftfx41(int n, double *a, int nw, double *w)
2693 {
2694     void cftf161(double *a, double *w);
2695     void cftf162(double *a, double *w);
2696     void cftf081(double *a, double *w);
2697     void cftf082(double *a, double *w);
2698
2699     if (n == 128) {
2700         cftf161(a, &w[nw - 8]);
2701         cftf162(&a[32], &w[nw - 32]);
2702         cftf161(&a[64], &w[nw - 8]);
2703         cftf161(&a[96], &w[nw - 8]);
2704     } else {
2705         cftf081(a, &w[nw - 8]);
2706         cftf082(&a[16], &w[nw - 8]);
2707         cftf081(&a[32], &w[nw - 8]);
2708         cftf081(&a[48], &w[nw - 8]);
2709     }
2710 }
2711
2712
2713 void cftf161(double *a, double *w)
2714 {
2715     double wn4r, wk1r, wk1i,
2716         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2717         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2718         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2719         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2720         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2721
2722     wn4r = w[1];
2723     wk1r = w[2];
2724     wk1i = w[3];
2725     x0r = a[0] + a[16];
2726     x0i = a[1] + a[17];
2727     x1r = a[0] - a[16];
2728     x1i = a[1] - a[17];
2729     x2r = a[8] + a[24];
2730     x2i = a[9] + a[25];
2731     x3r = a[8] - a[24];
2732     x3i = a[9] - a[25];
2733     y0r = x0r + x2r;
2734     y0i = x0i + x2i;
2735     y4r = x0r - x2r;
2736     y4i = x0i - x2i;
2737     y8r = x1r - x3i;
2738     y8i = x1i + x3r;
2739     y12r = x1r + x3i;
2740     y12i = x1i - x3r;
2741     x0r = a[2] + a[18];
2742     x0i = a[3] + a[19];
2743     x1r = a[2] - a[18];
2744     x1i = a[3] - a[19];
2745     x2r = a[10] + a[26];
2746     x2i = a[11] + a[27];
2747     x3r = a[10] - a[26];
2748     x3i = a[11] - a[27];
2749     y1r = x0r + x2r;
2750     y1i = x0i + x2i;
2751     y5r = x0r - x2r;
2752     y5i = x0i - x2i;
2753     x0r = x1r - x3i;
2754     x0i = x1i + x3r;
2755     y9r = wk1r * x0r - wk1i * x0i;
2756     y9i = wk1r * x0i + wk1i * x0r;
2757     x0r = x1r + x3i;
2758     x0i = x1i - x3r;
2759     y13r = wk1i * x0r - wk1r * x0i;
2760     y13i = wk1i * x0i + wk1r * x0r;
2761     x0r = a[4] + a[20];
2762     x0i = a[5] + a[21];
2763     x1r = a[4] - a[20];
2764     x1i = a[5] - a[21];
2765     x2r = a[12] + a[28];
2766     x2i = a[13] + a[29];
2767     x3r = a[12] - a[28];
2768     x3i = a[13] - a[29];
2769     y2r = x0r + x2r;
2770     y2i = x0i + x2i;
2771     y6r = x0r - x2r;
2772     y6i = x0i - x2i;
2773     x0r = x1r - x3i;
2774     x0i = x1i + x3r;
2775     y10r = wn4r * (x0r - x0i);
2776     y10i = wn4r * (x0i + x0r);
2777     x0r = x1r + x3i;
2778     x0i = x1i - x3r;
2779     y14r = wn4r * (x0r + x0i);
2780     y14i = wn4r * (x0i - x0r);
2781     x0r = a[6] + a[22];
2782     x0i = a[7] + a[23];
2783     x1r = a[6] - a[22];
2784     x1i = a[7] - a[23];
2785     x2r = a[14] + a[30];
2786     x2i = a[15] + a[31];
2787     x3r = a[14] - a[30];
2788     x3i = a[15] - a[31];
2789     y3r = x0r + x2r;
2790     y3i = x0i + x2i;
2791     y7r = x0r - x2r;
2792     y7i = x0i - x2i;
2793     x0r = x1r - x3i;
2794     x0i = x1i + x3r;
2795     y11r = wk1i * x0r - wk1r * x0i;
2796     y11i = wk1i * x0i + wk1r * x0r;
2797     x0r = x1r + x3i;
2798     x0i = x1i - x3r;
2799     y15r = wk1r * x0r - wk1i * x0i;
2800     y15i = wk1r * x0i + wk1i * x0r;
2801     x0r = y12r - y14r;
2802     x0i = y12i - y14i;
2803     x1r = y12r + y14r;
2804     x1i = y12i + y14i;
2805     x2r = y13r - y15r;
2806     x2i = y13i - y15i;
2807     x3r = y13r + y15r;
2808     x3i = y13i + y15i;
2809     a[24] = x0r + x2r;
2810     a[25] = x0i + x2i;
2811     a[26] = x0r - x2r;
2812     a[27] = x0i - x2i;
2813     a[28] = x1r - x3i;
2814     a[29] = x1i + x3r;
2815     a[30] = x1r + x3i;
2816     a[31] = x1i - x3r;
2817     x0r = y8r + y10r;
2818     x0i = y8i + y10i;
2819     x1r = y8r - y10r;
2820     x1i = y8i - y10i;
2821     x2r = y9r + y11r;
2822     x2i = y9i + y11i;
2823     x3r = y9r - y11r;
2824     x3i = y9i - y11i;
2825     a[16] = x0r + x2r;
2826     a[17] = x0i + x2i;
2827     a[18] = x0r - x2r;
2828     a[19] = x0i - x2i;
2829     a[20] = x1r - x3i;
2830     a[21] = x1i + x3r;
2831     a[22] = x1r + x3i;
2832     a[23] = x1i - x3r;
2833     x0r = y5r - y7i;
2834     x0i = y5i + y7r;
2835     x2r = wn4r * (x0r - x0i);
2836     x2i = wn4r * (x0i + x0r);
2837     x0r = y5r + y7i;
2838     x0i = y5i - y7r;
2839     x3r = wn4r * (x0r - x0i);
2840     x3i = wn4r * (x0i + x0r);
2841     x0r = y4r - y6i;
2842     x0i = y4i + y6r;
2843     x1r = y4r + y6i;
2844     x1i = y4i - y6r;
2845     a[8] = x0r + x2r;
2846     a[9] = x0i + x2i;
2847     a[10] = x0r - x2r;
2848     a[11] = x0i - x2i;
2849     a[12] = x1r - x3i;
2850     a[13] = x1i + x3r;
2851     a[14] = x1r + x3i;
2852     a[15] = x1i - x3r;
2853     x0r = y0r + y2r;
2854     x0i = y0i + y2i;
2855     x1r = y0r - y2r;
2856     x1i = y0i - y2i;
2857     x2r = y1r + y3r;
2858     x2i = y1i + y3i;
2859     x3r = y1r - y3r;
2860     x3i = y1i - y3i;
2861     a[0] = x0r + x2r;
2862     a[1] = x0i + x2i;
2863     a[2] = x0r - x2r;
2864     a[3] = x0i - x2i;
2865     a[4] = x1r - x3i;
2866     a[5] = x1i + x3r;
2867     a[6] = x1r + x3i;
2868     a[7] = x1i - x3r;
2869 }
2870
2871
2872 void cftf162(double *a, double *w)
2873 {
2874     double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
2875         x0r, x0i, x1r, x1i, x2r, x2i,
2876         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2877         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2878         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2879         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2880
2881     wn4r = w[1];
2882     wk1r = w[4];
2883     wk1i = w[5];
2884     wk3r = w[6];
2885     wk3i = -w[7];
2886     wk2r = w[8];
2887     wk2i = w[9];
2888     x1r = a[0] - a[17];
2889     x1i = a[1] + a[16];
2890     x0r = a[8] - a[25];
2891     x0i = a[9] + a[24];
2892     x2r = wn4r * (x0r - x0i);
2893     x2i = wn4r * (x0i + x0r);
2894     y0r = x1r + x2r;
2895     y0i = x1i + x2i;
2896     y4r = x1r - x2r;
2897     y4i = x1i - x2i;
2898     x1r = a[0] + a[17];
2899     x1i = a[1] - a[16];
2900     x0r = a[8] + a[25];
2901     x0i = a[9] - a[24];
2902     x2r = wn4r * (x0r - x0i);
2903     x2i = wn4r * (x0i + x0r);
2904     y8r = x1r - x2i;
2905     y8i = x1i + x2r;
2906     y12r = x1r + x2i;
2907     y12i = x1i - x2r;
2908     x0r = a[2] - a[19];
2909     x0i = a[3] + a[18];
2910     x1r = wk1r * x0r - wk1i * x0i;
2911     x1i = wk1r * x0i + wk1i * x0r;
2912     x0r = a[10] - a[27];
2913     x0i = a[11] + a[26];
2914     x2r = wk3i * x0r - wk3r * x0i;
2915     x2i = wk3i * x0i + wk3r * x0r;
2916     y1r = x1r + x2r;
2917     y1i = x1i + x2i;
2918     y5r = x1r - x2r;
2919     y5i = x1i - x2i;
2920     x0r = a[2] + a[19];
2921     x0i = a[3] - a[18];
2922     x1r = wk3r * x0r - wk3i * x0i;
2923     x1i = wk3r * x0i + wk3i * x0r;
2924     x0r = a[10] + a[27];
2925     x0i = a[11] - a[26];
2926     x2r = wk1r * x0r + wk1i * x0i;
2927     x2i = wk1r * x0i - wk1i * x0r;
2928     y9r = x1r - x2r;
2929     y9i = x1i - x2i;
2930     y13r = x1r + x2r;
2931     y13i = x1i + x2i;
2932     x0r = a[4] - a[21];
2933     x0i = a[5] + a[20];
2934     x1r = wk2r * x0r - wk2i * x0i;
2935     x1i = wk2r * x0i + wk2i * x0r;
2936     x0r = a[12] - a[29];
2937     x0i = a[13] + a[28];
2938     x2r = wk2i * x0r - wk2r * x0i;
2939     x2i = wk2i * x0i + wk2r * x0r;
2940     y2r = x1r + x2r;
2941     y2i = x1i + x2i;
2942     y6r = x1r - x2r;
2943     y6i = x1i - x2i;
2944     x0r = a[4] + a[21];
2945     x0i = a[5] - a[20];
2946     x1r = wk2i * x0r - wk2r * x0i;
2947     x1i = wk2i * x0i + wk2r * x0r;
2948     x0r = a[12] + a[29];
2949     x0i = a[13] - a[28];
2950     x2r = wk2r * x0r - wk2i * x0i;
2951     x2i = wk2r * x0i + wk2i * x0r;
2952     y10r = x1r - x2r;
2953     y10i = x1i - x2i;
2954     y14r = x1r + x2r;
2955     y14i = x1i + x2i;
2956     x0r = a[6] - a[23];
2957     x0i = a[7] + a[22];
2958     x1r = wk3r * x0r - wk3i * x0i;
2959     x1i = wk3r * x0i + wk3i * x0r;
2960     x0r = a[14] - a[31];
2961     x0i = a[15] + a[30];
2962     x2r = wk1i * x0r - wk1r * x0i;
2963     x2i = wk1i * x0i + wk1r * x0r;
2964     y3r = x1r + x2r;
2965     y3i = x1i + x2i;
2966     y7r = x1r - x2r;
2967     y7i = x1i - x2i;
2968     x0r = a[6] + a[23];
2969     x0i = a[7] - a[22];
2970     x1r = wk1i * x0r + wk1r * x0i;
2971     x1i = wk1i * x0i - wk1r * x0r;
2972     x0r = a[14] + a[31];
2973     x0i = a[15] - a[30];
2974     x2r = wk3i * x0r - wk3r * x0i;
2975     x2i = wk3i * x0i + wk3r * x0r;
2976     y11r = x1r + x2r;
2977     y11i = x1i + x2i;
2978     y15r = x1r - x2r;
2979     y15i = x1i - x2i;
2980     x1r = y0r + y2r;
2981     x1i = y0i + y2i;
2982     x2r = y1r + y3r;
2983     x2i = y1i + y3i;
2984     a[0] = x1r + x2r;
2985     a[1] = x1i + x2i;
2986     a[2] = x1r - x2r;
2987     a[3] = x1i - x2i;
2988     x1r = y0r - y2r;
2989     x1i = y0i - y2i;
2990     x2r = y1r - y3r;
2991     x2i = y1i - y3i;
2992     a[4] = x1r - x2i;
2993     a[5] = x1i + x2r;
2994     a[6] = x1r + x2i;
2995     a[7] = x1i - x2r;
2996     x1r = y4r - y6i;
2997     x1i = y4i + y6r;
2998     x0r = y5r - y7i;
2999     x0i = y5i + y7r;
3000     x2r = wn4r * (x0r - x0i);
3001     x2i = wn4r * (x0i + x0r);
3002     a[8] = x1r + x2r;
3003     a[9] = x1i + x2i;
3004     a[10] = x1r - x2r;
3005     a[11] = x1i - x2i;
3006     x1r = y4r + y6i;
3007     x1i = y4i - y6r;
3008     x0r = y5r + y7i;
3009     x0i = y5i - y7r;
3010     x2r = wn4r * (x0r - x0i);
3011     x2i = wn4r * (x0i + x0r);
3012     a[12] = x1r - x2i;
3013     a[13] = x1i + x2r;
3014     a[14] = x1r + x2i;
3015     a[15] = x1i - x2r;
3016     x1r = y8r + y10r;
3017     x1i = y8i + y10i;
3018     x2r = y9r - y11r;
3019     x2i = y9i - y11i;
3020     a[16] = x1r + x2r;
3021     a[17] = x1i + x2i;
3022     a[18] = x1r - x2r;
3023     a[19] = x1i - x2i;
3024     x1r = y8r - y10r;
3025     x1i = y8i - y10i;
3026     x2r = y9r + y11r;
3027     x2i = y9i + y11i;
3028     a[20] = x1r - x2i;
3029     a[21] = x1i + x2r;
3030     a[22] = x1r + x2i;
3031     a[23] = x1i - x2r;
3032     x1r = y12r - y14i;
3033     x1i = y12i + y14r;
3034     x0r = y13r + y15i;
3035     x0i = y13i - y15r;
3036     x2r = wn4r * (x0r - x0i);
3037     x2i = wn4r * (x0i + x0r);
3038     a[24] = x1r + x2r;
3039     a[25] = x1i + x2i;
3040     a[26] = x1r - x2r;
3041     a[27] = x1i - x2i;
3042     x1r = y12r + y14i;
3043     x1i = y12i - y14r;
3044     x0r = y13r - y15i;
3045     x0i = y13i + y15r;
3046     x2r = wn4r * (x0r - x0i);
3047     x2i = wn4r * (x0i + x0r);
3048     a[28] = x1r - x2i;
3049     a[29] = x1i + x2r;
3050     a[30] = x1r + x2i;
3051     a[31] = x1i - x2r;
3052 }
3053
3054
3055 void cftf081(double *a, double *w)
3056 {
3057     double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
3058         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3059         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3060
3061     wn4r = w[1];
3062     x0r = a[0] + a[8];
3063     x0i = a[1] + a[9];
3064     x1r = a[0] - a[8];
3065     x1i = a[1] - a[9];
3066     x2r = a[4] + a[12];
3067     x2i = a[5] + a[13];
3068     x3r = a[4] - a[12];
3069     x3i = a[5] - a[13];
3070     y0r = x0r + x2r;
3071     y0i = x0i + x2i;
3072     y2r = x0r - x2r;
3073     y2i = x0i - x2i;
3074     y1r = x1r - x3i;
3075     y1i = x1i + x3r;
3076     y3r = x1r + x3i;
3077     y3i = x1i - x3r;
3078     x0r = a[2] + a[10];
3079     x0i = a[3] + a[11];
3080     x1r = a[2] - a[10];
3081     x1i = a[3] - a[11];
3082     x2r = a[6] + a[14];
3083     x2i = a[7] + a[15];
3084     x3r = a[6] - a[14];
3085     x3i = a[7] - a[15];
3086     y4r = x0r + x2r;
3087     y4i = x0i + x2i;
3088     y6r = x0r - x2r;
3089     y6i = x0i - x2i;
3090     x0r = x1r - x3i;
3091     x0i = x1i + x3r;
3092     x2r = x1r + x3i;
3093     x2i = x1i - x3r;
3094     y5r = wn4r * (x0r - x0i);
3095     y5i = wn4r * (x0r + x0i);
3096     y7r = wn4r * (x2r - x2i);
3097     y7i = wn4r * (x2r + x2i);
3098     a[8] = y1r + y5r;
3099     a[9] = y1i + y5i;
3100     a[10] = y1r - y5r;
3101     a[11] = y1i - y5i;
3102     a[12] = y3r - y7i;
3103     a[13] = y3i + y7r;
3104     a[14] = y3r + y7i;
3105     a[15] = y3i - y7r;
3106     a[0] = y0r + y4r;
3107     a[1] = y0i + y4i;
3108     a[2] = y0r - y4r;
3109     a[3] = y0i - y4i;
3110     a[4] = y2r - y6i;
3111     a[5] = y2i + y6r;
3112     a[6] = y2r + y6i;
3113     a[7] = y2i - y6r;
3114 }
3115
3116
3117 void cftf082(double *a, double *w)
3118 {
3119     double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
3120         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3121         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3122
3123     wn4r = w[1];
3124     wk1r = w[2];
3125     wk1i = w[3];
3126     y0r = a[0] - a[9];
3127     y0i = a[1] + a[8];
3128     y1r = a[0] + a[9];
3129     y1i = a[1] - a[8];
3130     x0r = a[4] - a[13];
3131     x0i = a[5] + a[12];
3132     y2r = wn4r * (x0r - x0i);
3133     y2i = wn4r * (x0i + x0r);
3134     x0r = a[4] + a[13];
3135     x0i = a[5] - a[12];
3136     y3r = wn4r * (x0r - x0i);
3137     y3i = wn4r * (x0i + x0r);
3138     x0r = a[2] - a[11];
3139     x0i = a[3] + a[10];
3140     y4r = wk1r * x0r - wk1i * x0i;
3141     y4i = wk1r * x0i + wk1i * x0r;
3142     x0r = a[2] + a[11];
3143     x0i = a[3] - a[10];
3144     y5r = wk1i * x0r - wk1r * x0i;
3145     y5i = wk1i * x0i + wk1r * x0r;
3146     x0r = a[6] - a[15];
3147     x0i = a[7] + a[14];
3148     y6r = wk1i * x0r - wk1r * x0i;
3149     y6i = wk1i * x0i + wk1r * x0r;
3150     x0r = a[6] + a[15];
3151     x0i = a[7] - a[14];
3152     y7r = wk1r * x0r - wk1i * x0i;
3153     y7i = wk1r * x0i + wk1i * x0r;
3154     x0r = y0r + y2r;
3155     x0i = y0i + y2i;
3156     x1r = y4r + y6r;
3157     x1i = y4i + y6i;
3158     a[0] = x0r + x1r;
3159     a[1] = x0i + x1i;
3160     a[2] = x0r - x1r;
3161     a[3] = x0i - x1i;
3162     x0r = y0r - y2r;
3163     x0i = y0i - y2i;
3164     x1r = y4r - y6r;
3165     x1i = y4i - y6i;
3166     a[4] = x0r - x1i;
3167     a[5] = x0i + x1r;
3168     a[6] = x0r + x1i;
3169     a[7] = x0i - x1r;
3170     x0r = y1r - y3i;
3171     x0i = y1i + y3r;
3172     x1r = y5r - y7r;
3173     x1i = y5i - y7i;
3174     a[8] = x0r + x1r;
3175     a[9] = x0i + x1i;
3176     a[10] = x0r - x1r;
3177     a[11] = x0i - x1i;
3178     x0r = y1r + y3i;
3179     x0i = y1i - y3r;
3180     x1r = y5r + y7r;
3181     x1i = y5i + y7i;
3182     a[12] = x0r - x1i;
3183     a[13] = x0i + x1r;
3184     a[14] = x0r + x1i;
3185     a[15] = x0i - x1r;
3186 }
3187
3188
3189 void cftf040(double *a)
3190 {
3191     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3192
3193     x0r = a[0] + a[4];
3194     x0i = a[1] + a[5];
3195     x1r = a[0] - a[4];
3196     x1i = a[1] - a[5];
3197     x2r = a[2] + a[6];
3198     x2i = a[3] + a[7];
3199     x3r = a[2] - a[6];
3200     x3i = a[3] - a[7];
3201     a[0] = x0r + x2r;
3202     a[1] = x0i + x2i;
3203     a[2] = x1r - x3i;
3204     a[3] = x1i + x3r;
3205     a[4] = x0r - x2r;
3206     a[5] = x0i - x2i;
3207     a[6] = x1r + x3i;
3208     a[7] = x1i - x3r;
3209 }
3210
3211
3212 void cftb040(double *a)
3213 {
3214     double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3215
3216     x0r = a[0] + a[4];
3217     x0i = a[1] + a[5];
3218     x1r = a[0] - a[4];
3219     x1i = a[1] - a[5];
3220     x2r = a[2] + a[6];
3221     x2i = a[3] + a[7];
3222     x3r = a[2] - a[6];
3223     x3i = a[3] - a[7];
3224     a[0] = x0r + x2r;
3225     a[1] = x0i + x2i;
3226     a[2] = x1r + x3i;
3227     a[3] = x1i - x3r;
3228     a[4] = x0r - x2r;
3229     a[5] = x0i - x2i;
3230     a[6] = x1r - x3i;
3231     a[7] = x1i + x3r;
3232 }
3233
3234
3235 void cftx020(double *a)
3236 {
3237     double x0r, x0i;
3238
3239     x0r = a[0] - a[2];
3240     x0i = a[1] - a[3];
3241     a[0] += a[2];
3242     a[1] += a[3];
3243     a[2] = x0r;
3244     a[3] = x0i;
3245 }
3246
3247
3248 void rftfsub(int n, double *a, int nc, double *c)
3249 {
3250     int j, k, kk, ks, m;
3251     double wkr, wki, xr, xi, yr, yi;
3252
3253     m = n >> 1;
3254     ks = 2 * nc / m;
3255     kk = 0;
3256     for (j = 2; j < m; j += 2) {
3257         k = n - j;
3258         kk += ks;
3259         wkr = 0.5 - c[nc - kk];
3260         wki = c[kk];
3261         xr = a[j] - a[k];
3262         xi = a[j + 1] + a[k + 1];
3263         yr = wkr * xr - wki * xi;
3264         yi = wkr * xi + wki * xr;
3265         a[j] -= yr;
3266         a[j + 1] -= yi;
3267         a[k] += yr;
3268         a[k + 1] -= yi;
3269     }
3270 }
3271
3272
3273 void rftbsub(int n, double *a, int nc, double *c)
3274 {
3275     int j, k, kk, ks, m;
3276     double wkr, wki, xr, xi, yr, yi;
3277
3278     m = n >> 1;
3279     ks = 2 * nc / m;
3280     kk = 0;
3281     for (j = 2; j < m; j += 2) {
3282         k = n - j;
3283         kk += ks;
3284         wkr = 0.5 - c[nc - kk];
3285         wki = c[kk];
3286         xr = a[j] - a[k];
3287         xi = a[j + 1] + a[k + 1];
3288         yr = wkr * xr + wki * xi;
3289         yi = wkr * xi - wki * xr;
3290         a[j] -= yr;
3291         a[j + 1] -= yi;
3292         a[k] += yr;
3293         a[k + 1] -= yi;
3294     }
3295 }
3296
3297
3298 void dctsub(int n, double *a, int nc, double *c)
3299 {
3300     int j, k, kk, ks, m;
3301     double wkr, wki, xr;
3302
3303     m = n >> 1;
3304     ks = nc / n;
3305     kk = 0;
3306     for (j = 1; j < m; j++) {
3307         k = n - j;
3308         kk += ks;
3309         wkr = c[kk] - c[nc - kk];
3310         wki = c[kk] + c[nc - kk];
3311         xr = wki * a[j] - wkr * a[k];
3312         a[j] = wkr * a[j] + wki * a[k];
3313         a[k] = xr;
3314     }
3315     a[m] *= c[0];
3316 }
3317
3318
3319 void dstsub(int n, double *a, int nc, double *c)
3320 {
3321     int j, k, kk, ks, m;
3322     double wkr, wki, xr;
3323
3324     m = n >> 1;
3325     ks = nc / n;
3326     kk = 0;
3327     for (j = 1; j < m; j++) {
3328         k = n - j;
3329         kk += ks;
3330         wkr = c[kk] - c[nc - kk];
3331         wki = c[kk] + c[nc - kk];
3332         xr = wki * a[k] - wkr * a[j];
3333         a[k] = wkr * a[k] + wki * a[j];
3334         a[j] = xr;
3335     }
3336     a[m] *= c[0];
3337 }
3338