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

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