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

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