src/mdlib/fftpack.c

   1 /*
   2
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2012, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Groningen Machine for Chemical Simulation
  33
  34 ************************************************************
  35
  36 Copy of fftpack from Numpy with very minor modifications:
  37   - usage of fftpack.h (replacement for Treal define)
  38   - [cr]fft[ifb]1 non-static
  39   - Added Copyright headers
  40   - Added fftpack_ prefix
  41
  42 Original version is from Numpy 1.6
  43
  44 ************************************************************
  45
  46 Copyright (c) 2005-2011, NumPy Developers.
  47 All rights reserved.
  48
  49 Redistribution and use in source and binary forms, with or without
  50 modification, are permitted provided that the following conditions are
  51 met:
  52
  53     * Redistributions of source code must retain the above copyright
  54        notice, this list of conditions and the following disclaimer.
  55
  56     * Redistributions in binary form must reproduce the above
  57        copyright notice, this list of conditions and the following
  58        disclaimer in the documentation and/or other materials provided
  59        with the distribution.
  60
  61     * Neither the name of the NumPy Developers nor the names of any
  62        contributors may be used to endorse or promote products derived
  63        from this software without specific prior written permission.
  64
  65 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  66 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  67 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  68 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  69 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  70 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  71 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  72 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  73 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  74 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  75 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  76
  77 fftpack.c : A set of FFT routines in C.
  78 Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
  79
  80 */
  81
  82 /* isign is +1 for backward and -1 for forward transforms */
  83
  84 #include <math.h>
  85 #include <stdio.h>
  86
  87 #include "fftpack.h"
  88
  89 #define ref(u,a) u[a]
  90
  91 #define MAXFAC 13    /* maximum number of factors in factorization of n */
  92 #define NSPECIAL 4   /* number of factors for which we have special-case routines */
  93
  94 #ifdef __cplusplus
  95 extern "C" {
  96 #endif
  97
  98
  99 /* ----------------------------------------------------------------------
 100    passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
 101 ---------------------------------------------------------------------- */
 102
 103 static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign)
 104   /* isign==+1 for backward transform */
 105   {
 106     int i, k, ah, ac;
 107     Treal ti2, tr2;
 108     if (ido <= 2) {
 109       for (k=0; k<l1; k++) {
 110         ah = k*ido;
 111         ac = 2*k*ido;
 112         ch[ah]              = ref(cc,ac) + ref(cc,ac + ido);
 113         ch[ah + ido*l1]     = ref(cc,ac) - ref(cc,ac + ido);
 114         ch[ah+1]            = ref(cc,ac+1) + ref(cc,ac + ido + 1);
 115         ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1);
 116       }
 117     } else {
 118       for (k=0; k<l1; k++) {
 119         for (i=0; i<ido-1; i+=2) {
 120           ah = i + k*ido;
 121           ac = i + 2*k*ido;
 122           ch[ah]   = ref(cc,ac) + ref(cc,ac + ido);
 123           tr2      = ref(cc,ac) - ref(cc,ac + ido);
 124           ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido);
 125           ti2      = ref(cc,ac+1) - ref(cc,ac + 1 + ido);
 126           ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
 127           ch[ah+l1*ido]   = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
 128         }
 129       }
 130     }
 131   } /* passf2 */
 132
 133
 134 static void passf3(int ido, int l1, const Treal cc[], Treal ch[],
 135       const Treal wa1[], const Treal wa2[], int isign)
 136   /* isign==+1 for backward transform */
 137   {
 138     static const Treal taur = -0.5;
 139     static const Treal taui = 0.866025403784439;
 140     int i, k, ac, ah;
 141     Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
 142     if (ido == 2) {
 143       for (k=1; k<=l1; k++) {
 144         ac = (3*k - 2)*ido;
 145         tr2 = ref(cc,ac) + ref(cc,ac + ido);
 146         cr2 = ref(cc,ac - ido) + taur*tr2;
 147         ah = (k - 1)*ido;
 148         ch[ah] = ref(cc,ac - ido) + tr2;
 149
 150         ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
 151         ci2 = ref(cc,ac - ido + 1) + taur*ti2;
 152         ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
 153
 154         cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
 155         ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
 156         ch[ah + l1*ido] = cr2 - ci3;
 157         ch[ah + 2*l1*ido] = cr2 + ci3;
 158         ch[ah + l1*ido + 1] = ci2 + cr3;
 159         ch[ah + 2*l1*ido + 1] = ci2 - cr3;
 160       }
 161     } else {
 162       for (k=1; k<=l1; k++) {
 163         for (i=0; i<ido-1; i+=2) {
 164           ac = i + (3*k - 2)*ido;
 165           tr2 = ref(cc,ac) + ref(cc,ac + ido);
 166           cr2 = ref(cc,ac - ido) + taur*tr2;
 167           ah = i + (k-1)*ido;
 168           ch[ah] = ref(cc,ac - ido) + tr2;
 169           ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1);
 170           ci2 = ref(cc,ac - ido + 1) + taur*ti2;
 171           ch[ah + 1] = ref(cc,ac - ido + 1) + ti2;
 172           cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido));
 173           ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1));
 174           dr2 = cr2 - ci3;
 175           dr3 = cr2 + ci3;
 176           di2 = ci2 + cr3;
 177           di3 = ci2 - cr3;
 178           ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
 179           ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
 180           ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
 181           ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
 182         }
 183       }
 184     }
 185   } /* passf3 */
 186
 187
 188 static void passf4(int ido, int l1, const Treal cc[], Treal ch[],
 189       const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign)
 190   /* isign == -1 for forward transform and +1 for backward transform */
 191   {
 192     int i, k, ac, ah;
 193     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
 194     if (ido == 2) {
 195       for (k=0; k<l1; k++) {
 196         ac = 4*k*ido + 1;
 197         ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
 198         ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
 199         tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
 200         ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
 201         tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
 202         tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
 203         ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
 204         tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
 205         ah = k*ido;
 206         ch[ah] = tr2 + tr3;
 207         ch[ah + 2*l1*ido] = tr2 - tr3;
 208         ch[ah + 1] = ti2 + ti3;
 209         ch[ah + 2*l1*ido + 1] = ti2 - ti3;
 210         ch[ah + l1*ido] = tr1 + isign*tr4;
 211         ch[ah + 3*l1*ido] = tr1 - isign*tr4;
 212         ch[ah + l1*ido + 1] = ti1 + isign*ti4;
 213         ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4;
 214       }
 215     } else {
 216       for (k=0; k<l1; k++) {
 217         for (i=0; i<ido-1; i+=2) {
 218           ac = i + 1 + 4*k*ido;
 219           ti1 = ref(cc,ac) - ref(cc,ac + 2*ido);
 220           ti2 = ref(cc,ac) + ref(cc,ac + 2*ido);
 221           ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido);
 222           tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido);
 223           tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1);
 224           tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1);
 225           ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1);
 226           tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1);
 227           ah = i + k*ido;
 228           ch[ah] = tr2 + tr3;
 229           cr3 = tr2 - tr3;
 230           ch[ah + 1] = ti2 + ti3;
 231           ci3 = ti2 - ti3;
 232           cr2 = tr1 + isign*tr4;
 233           cr4 = tr1 - isign*tr4;
 234           ci2 = ti1 + isign*ti4;
 235           ci4 = ti1 - isign*ti4;
 236           ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
 237           ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
 238           ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
 239           ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
 240           ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
 241           ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
 242         }
 243       }
 244     }
 245   } /* passf4 */
 246
 247
 248 static void passf5(int ido, int l1, const Treal cc[], Treal ch[],
 249       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign)
 250   /* isign == -1 for forward transform and +1 for backward transform */
 251   {
 252     static const Treal tr11 = 0.309016994374947;
 253     static const Treal ti11 = 0.951056516295154;
 254     static const Treal tr12 = -0.809016994374947;
 255     static const Treal ti12 = 0.587785252292473;
 256     int i, k, ac, ah;
 257     Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
 258         ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
 259     if (ido == 2) {
 260       for (k = 1; k <= l1; ++k) {
 261         ac = (5*k - 4)*ido + 1;
 262         ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
 263         ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
 264         ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
 265         ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
 266         tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
 267         tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
 268         tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
 269         tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
 270         ah = (k - 1)*ido;
 271         ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
 272         ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
 273         cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
 274         ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
 275         cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
 276         ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
 277         cr5 = isign*(ti11*tr5 + ti12*tr4);
 278         ci5 = isign*(ti11*ti5 + ti12*ti4);
 279         cr4 = isign*(ti12*tr5 - ti11*tr4);
 280         ci4 = isign*(ti12*ti5 - ti11*ti4);
 281         ch[ah + l1*ido] = cr2 - ci5;
 282         ch[ah + 4*l1*ido] = cr2 + ci5;
 283         ch[ah + l1*ido + 1] = ci2 + cr5;
 284         ch[ah + 2*l1*ido + 1] = ci3 + cr4;
 285         ch[ah + 2*l1*ido] = cr3 - ci4;
 286         ch[ah + 3*l1*ido] = cr3 + ci4;
 287         ch[ah + 3*l1*ido + 1] = ci3 - cr4;
 288         ch[ah + 4*l1*ido + 1] = ci2 - cr5;
 289       }
 290     } else {
 291       for (k=1; k<=l1; k++) {
 292         for (i=0; i<ido-1; i+=2) {
 293           ac = i + 1 + (k*5 - 4)*ido;
 294           ti5 = ref(cc,ac) - ref(cc,ac + 3*ido);
 295           ti2 = ref(cc,ac) + ref(cc,ac + 3*ido);
 296           ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido);
 297           ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido);
 298           tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1);
 299           tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1);
 300           tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1);
 301           tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1);
 302           ah = i + (k - 1)*ido;
 303           ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3;
 304           ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3;
 305           cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
 306
 307           ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3;
 308           cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
 309
 310           ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3;
 311           cr5 = isign*(ti11*tr5 + ti12*tr4);
 312           ci5 = isign*(ti11*ti5 + ti12*ti4);
 313           cr4 = isign*(ti12*tr5 - ti11*tr4);
 314           ci4 = isign*(ti12*ti5 - ti11*ti4);
 315           dr3 = cr3 - ci4;
 316           dr4 = cr3 + ci4;
 317           di3 = ci3 + cr4;
 318           di4 = ci3 - cr4;
 319           dr5 = cr2 + ci5;
 320           dr2 = cr2 - ci5;
 321           di5 = ci2 - cr5;
 322           di2 = ci2 + cr5;
 323           ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2;
 324           ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2;
 325           ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3;
 326           ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3;
 327           ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4;
 328           ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4;
 329           ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5;
 330           ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5;
 331         }
 332       }
 333     }
 334   } /* passf5 */
 335
 336
 337 static void passf(int *nac, int ido, int ip, int l1, int idl1,
 338       Treal cc[], Treal ch[],
 339       const Treal wa[], int isign)
 340   /* isign is -1 for forward transform and +1 for backward transform */
 341   {
 342     int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp;
 343     Treal wai, war;
 344
 345     idot = ido / 2;
 346     /* nt = ip*idl1;*/
 347     ipph = (ip + 1) / 2;
 348     idp = ip*ido;
 349     if (ido >= l1) {
 350       for (j=1; j<ipph; j++) {
 351         jc = ip - j;
 352         for (k=0; k<l1; k++) {
 353           for (i=0; i<ido; i++) {
 354             ch[i + (k + j*l1)*ido] =
 355                 ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido);
 356             ch[i + (k + jc*l1)*ido] =
 357                 ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido);
 358           }
 359         }
 360       }
 361       for (k=0; k<l1; k++)
 362         for (i=0; i<ido; i++)
 363           ch[i + k*ido] = ref(cc,i + k*ip*ido);
 364     } else {
 365       for (j=1; j<ipph; j++) {
 366         jc = ip - j;
 367         for (i=0; i<ido; i++) {
 368           for (k=0; k<l1; k++) {
 369             ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*
 370                 ip)*ido);
 371             ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*
 372                 ip)*ido);
 373           }
 374         }
 375       }
 376       for (i=0; i<ido; i++)
 377         for (k=0; k<l1; k++)
 378           ch[i + k*ido] = ref(cc,i + k*ip*ido);
 379     }
 380
 381     idl = 2 - ido;
 382     inc = 0;
 383     for (l=1; l<ipph; l++) {
 384       lc = ip - l;
 385       idl += ido;
 386       for (ik=0; ik<idl1; ik++) {
 387         cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1];
 388         cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1];
 389       }
 390       idlj = idl;
 391       inc += ido;
 392       for (j=2; j<ipph; j++) {
 393         jc = ip - j;
 394         idlj += inc;
 395         if (idlj > idp) idlj -= idp;
 396         war = wa[idlj - 2];
 397         wai = wa[idlj-1];
 398         for (ik=0; ik<idl1; ik++) {
 399           cc[ik + l*idl1] += war*ch[ik + j*idl1];
 400           cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1];
 401         }
 402       }
 403     }
 404     for (j=1; j<ipph; j++)
 405       for (ik=0; ik<idl1; ik++)
 406         ch[ik] += ch[ik + j*idl1];
 407     for (j=1; j<ipph; j++) {
 408       jc = ip - j;
 409       for (ik=1; ik<idl1; ik+=2) {
 410         ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1];
 411         ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1];
 412         ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1];
 413         ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1];
 414       }
 415     }
 416     *nac = 1;
 417     if (ido == 2) return;
 418     *nac = 0;
 419     for (ik=0; ik<idl1; ik++)
 420       cc[ik] = ch[ik];
 421     for (j=1; j<ip; j++) {
 422       for (k=0; k<l1; k++) {
 423         cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0];
 424         cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1];
 425       }
 426     }
 427     if (idot <= l1) {
 428       idij = 0;
 429       for (j=1; j<ip; j++) {
 430         idij += 2;
 431         for (i=3; i<ido; i+=2) {
 432           idij += 2;
 433           for (k=0; k<l1; k++) {
 434             cc[i - 1 + (k + j*l1)*ido] =
 435                 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
 436                 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
 437             cc[i + (k + j*l1)*ido] =
 438                 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
 439                 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
 440           }
 441         }
 442       }
 443     } else {
 444       idj = 2 - ido;
 445       for (j=1; j<ip; j++) {
 446         idj += ido;
 447         for (k = 0; k < l1; k++) {
 448           idij = idj;
 449           for (i=3; i<ido; i+=2) {
 450             idij += 2;
 451             cc[i - 1 + (k + j*l1)*ido] =
 452                 wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] -
 453                 isign*wa[idij-1]*ch[i + (k + j*l1)*ido];
 454             cc[i + (k + j*l1)*ido] =
 455                 wa[idij - 2]*ch[i + (k + j*l1)*ido] +
 456                 isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido];
 457           }
 458         }
 459       }
 460     }
 461   } /* passf */
 462
 463
 464   /* ----------------------------------------------------------------------
 465 radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
 466 Treal FFT passes fwd and bwd.
 467 ---------------------------------------------------------------------- */
 468
 469 static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
 470   {
 471     int i, k, ic;
 472     Treal ti2, tr2;
 473     for (k=0; k<l1; k++) {
 474       ch[2*k*ido] =
 475           ref(cc,k*ido) + ref(cc,(k + l1)*ido);
 476       ch[(2*k+1)*ido + ido-1] =
 477           ref(cc,k*ido) - ref(cc,(k + l1)*ido);
 478     }
 479     if (ido < 2) return;
 480     if (ido != 2) {
 481       for (k=0; k<l1; k++) {
 482         for (i=2; i<ido; i+=2) {
 483           ic = ido - i;
 484           tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
 485           ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
 486           ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2;
 487           ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido);
 488           ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2;
 489           ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2;
 490         }
 491       }
 492       if (ido % 2 == 1) return;
 493     }
 494     for (k=0; k<l1; k++) {
 495       ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido);
 496       ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido);
 497     }
 498   } /* radf2 */
 499
 500
 501 static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[])
 502   {
 503     int i, k, ic;
 504     Treal ti2, tr2;
 505     for (k=0; k<l1; k++) {
 506       ch[k*ido] =
 507           ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
 508       ch[(k + l1)*ido] =
 509           ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
 510     }
 511     if (ido < 2) return;
 512     if (ido != 2) {
 513       for (k = 0; k < l1; ++k) {
 514         for (i = 2; i < ido; i += 2) {
 515           ic = ido - i;
 516           ch[i-1 + k*ido] =
 517               ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
 518           tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
 519           ch[i + k*ido] =
 520               ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
 521           ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
 522           ch[i-1 + (k + l1)*ido] =
 523               wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
 524           ch[i + (k + l1)*ido] =
 525               wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
 526         }
 527       }
 528       if (ido % 2 == 1) return;
 529     }
 530     for (k = 0; k < l1; k++) {
 531       ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido);
 532       ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido);
 533     }
 534   } /* radb2 */
 535
 536
 537 static void radf3(int ido, int l1, const Treal cc[], Treal ch[],
 538       const Treal wa1[], const Treal wa2[])
 539   {
 540     static const Treal taur = -0.5;
 541     static const Treal taui = 0.866025403784439;
 542     int i, k, ic;
 543     Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
 544     for (k=0; k<l1; k++) {
 545       cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
 546       ch[3*k*ido] = ref(cc,k*ido) + cr2;
 547       ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
 548       ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2;
 549     }
 550     if (ido == 1) return;
 551     for (k=0; k<l1; k++) {
 552       for (i=2; i<ido; i+=2) {
 553         ic = ido - i;
 554         dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) +
 555             wa1[i - 1]*ref(cc,i + (k + l1)*ido);
 556         di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
 557         dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
 558         di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
 559         cr2 = dr2 + dr3;
 560         ci2 = di2 + di3;
 561         ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2;
 562         ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2;
 563         tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
 564         ti2 = ref(cc,i + k*ido) + taur*ci2;
 565         tr3 = taui*(di2 - di3);
 566         ti3 = taui*(dr3 - dr2);
 567         ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3;
 568         ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3;
 569         ch[i + (3*k + 2)*ido] = ti2 + ti3;
 570         ch[ic + (3*k + 1)*ido] = ti3 - ti2;
 571       }
 572     }
 573   } /* radf3 */
 574
 575
 576 static void radb3(int ido, int l1, const Treal cc[], Treal ch[],
 577       const Treal wa1[], const Treal wa2[])
 578   {
 579     static const Treal taur = -0.5;
 580     static const Treal taui = 0.866025403784439;
 581     int i, k, ic;
 582     Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
 583     for (k=0; k<l1; k++) {
 584       tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
 585       cr2 = ref(cc,3*k*ido) + taur*tr2;
 586       ch[k*ido] = ref(cc,3*k*ido) + tr2;
 587       ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
 588       ch[(k + l1)*ido] = cr2 - ci3;
 589       ch[(k + 2*l1)*ido] = cr2 + ci3;
 590     }
 591     if (ido == 1) return;
 592     for (k=0; k<l1; k++) {
 593       for (i=2; i<ido; i+=2) {
 594         ic = ido - i;
 595         tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
 596         cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
 597         ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2;
 598         ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
 599         ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
 600         ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2;
 601         cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
 602         ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
 603         dr2 = cr2 - ci3;
 604         dr3 = cr2 + ci3;
 605         di2 = ci2 + cr3;
 606         di3 = ci2 - cr3;
 607         ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
 608         ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
 609         ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
 610         ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
 611       }
 612     }
 613   } /* radb3 */
 614
 615
 616 static void radf4(int ido, int l1, const Treal cc[], Treal ch[],
 617       const Treal wa1[], const Treal wa2[], const Treal wa3[])
 618   {
 619     static const Treal hsqt2 = 0.7071067811865475;
 620     int i, k, ic;
 621     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
 622     for (k=0; k<l1; k++) {
 623       tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
 624       tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
 625       ch[4*k*ido] = tr1 + tr2;
 626       ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1;
 627       ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
 628       ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
 629     }
 630     if (ido < 2) return;
 631     if (ido != 2) {
 632       for (k=0; k<l1; k++) {
 633         for (i=2; i<ido; i += 2) {
 634           ic = ido - i;
 635           cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
 636           ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
 637           cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*
 638               ido);
 639           ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*
 640               ido);
 641           cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*
 642               ido);
 643           ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*
 644               ido);
 645           tr1 = cr2 + cr4;
 646           tr4 = cr4 - cr2;
 647           ti1 = ci2 + ci4;
 648           ti4 = ci2 - ci4;
 649           ti2 = ref(cc,i + k*ido) + ci3;
 650           ti3 = ref(cc,i + k*ido) - ci3;
 651           tr2 = ref(cc,i - 1 + k*ido) + cr3;
 652           tr3 = ref(cc,i - 1 + k*ido) - cr3;
 653           ch[i - 1 + 4*k*ido] = tr1 + tr2;
 654           ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1;
 655           ch[i + 4*k*ido] = ti1 + ti2;
 656           ch[ic + (4*k + 3)*ido] = ti1 - ti2;
 657           ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3;
 658           ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4;
 659           ch[i + (4*k + 2)*ido] = tr4 + ti3;
 660           ch[ic + (4*k + 1)*ido] = tr4 - ti3;
 661         }
 662       }
 663       if (ido % 2 == 1) return;
 664     }
 665     for (k=0; k<l1; k++) {
 666       ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
 667       tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
 668       ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido);
 669       ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1;
 670       ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
 671       ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
 672     }
 673   } /* radf4 */
 674
 675
 676 static void radb4(int ido, int l1, const Treal cc[], Treal ch[],
 677       const Treal wa1[], const Treal wa2[], const Treal wa3[])
 678   {
 679     static const Treal sqrt2 = 1.414213562373095;
 680     int i, k, ic;
 681     Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
 682     for (k = 0; k < l1; k++) {
 683       tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
 684       tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
 685       tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
 686       tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
 687       ch[k*ido] = tr2 + tr3;
 688       ch[(k + l1)*ido] = tr1 - tr4;
 689       ch[(k + 2*l1)*ido] = tr2 - tr3;
 690       ch[(k + 3*l1)*ido] = tr1 + tr4;
 691     }
 692     if (ido < 2) return;
 693     if (ido != 2) {
 694       for (k = 0; k < l1; ++k) {
 695         for (i = 2; i < ido; i += 2) {
 696           ic = ido - i;
 697           ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
 698           ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
 699           ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
 700           tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
 701           tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
 702           tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
 703           ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
 704           tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
 705           ch[i - 1 + k*ido] = tr2 + tr3;
 706           cr3 = tr2 - tr3;
 707           ch[i + k*ido] = ti2 + ti3;
 708           ci3 = ti2 - ti3;
 709           cr2 = tr1 - tr4;
 710           cr4 = tr1 + tr4;
 711           ci2 = ti1 + ti4;
 712           ci4 = ti1 - ti4;
 713           ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
 714           ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
 715           ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
 716           ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
 717           ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
 718           ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
 719         }
 720       }
 721       if (ido % 2 == 1) return;
 722     }
 723     for (k = 0; k < l1; k++) {
 724       ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
 725       ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
 726       tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
 727       tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
 728       ch[ido-1 + k*ido] = tr2 + tr2;
 729       ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1);
 730       ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2;
 731       ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1);
 732     }
 733   } /* radb4 */
 734
 735
 736 static void radf5(int ido, int l1, const Treal cc[], Treal ch[],
 737       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
 738   {
 739     static const Treal tr11 = 0.309016994374947;
 740     static const Treal ti11 = 0.951056516295154;
 741     static const Treal tr12 = -0.809016994374947;
 742     static const Treal ti12 = 0.587785252292473;
 743     int i, k, ic;
 744     Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
 745         cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
 746     for (k = 0; k < l1; k++) {
 747       cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
 748       ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
 749       cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
 750       ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
 751       ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3;
 752       ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
 753       ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4;
 754       ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
 755       ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4;
 756     }
 757     if (ido == 1) return;
 758     for (k = 0; k < l1; ++k) {
 759       for (i = 2; i < ido; i += 2) {
 760         ic = ido - i;
 761         dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
 762         di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
 763         dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
 764         di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
 765         dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
 766         di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
 767         dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
 768         di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
 769         cr2 = dr2 + dr5;
 770         ci5 = dr5 - dr2;
 771         cr5 = di2 - di5;
 772         ci2 = di2 + di5;
 773         cr3 = dr3 + dr4;
 774         ci4 = dr4 - dr3;
 775         cr4 = di3 - di4;
 776         ci3 = di3 + di4;
 777         ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3;
 778         ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3;
 779         tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
 780         ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
 781         tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
 782         ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
 783         tr5 = ti11*cr5 + ti12*cr4;
 784         ti5 = ti11*ci5 + ti12*ci4;
 785         tr4 = ti12*cr5 - ti11*cr4;
 786         ti4 = ti12*ci5 - ti11*ci4;
 787         ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5;
 788         ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5;
 789         ch[i + (5*k + 2)*ido] = ti2 + ti5;
 790         ch[ic + (5*k + 1)*ido] = ti5 - ti2;
 791         ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4;
 792         ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4;
 793         ch[i + (5*k + 4)*ido] = ti3 + ti4;
 794         ch[ic + (5*k + 3)*ido] = ti4 - ti3;
 795       }
 796     }
 797   } /* radf5 */
 798
 799
 800 static void radb5(int ido, int l1, const Treal cc[], Treal ch[],
 801       const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[])
 802   {
 803     static const Treal tr11 = 0.309016994374947;
 804     static const Treal ti11 = 0.951056516295154;
 805     static const Treal tr12 = -0.809016994374947;
 806     static const Treal ti12 = 0.587785252292473;
 807     int i, k, ic;
 808     Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
 809         ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
 810     for (k = 0; k < l1; k++) {
 811       ti5 = 2*ref(cc,(5*k + 2)*ido);
 812       ti4 = 2*ref(cc,(5*k + 4)*ido);
 813       tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
 814       tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
 815       ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3;
 816       cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
 817       cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
 818       ci5 = ti11*ti5 + ti12*ti4;
 819       ci4 = ti12*ti5 - ti11*ti4;
 820       ch[(k + l1)*ido] = cr2 - ci5;
 821       ch[(k + 2*l1)*ido] = cr3 - ci4;
 822       ch[(k + 3*l1)*ido] = cr3 + ci4;
 823       ch[(k + 4*l1)*ido] = cr2 + ci5;
 824     }
 825     if (ido == 1) return;
 826     for (k = 0; k < l1; ++k) {
 827       for (i = 2; i < ido; i += 2) {
 828         ic = ido - i;
 829         ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
 830         ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
 831         ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
 832         ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
 833         tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
 834         tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
 835         tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
 836         tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
 837         ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
 838         ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3;
 839         cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;
 840
 841         ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
 842         cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;
 843
 844         ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3;
 845         cr5 = ti11*tr5 + ti12*tr4;
 846         ci5 = ti11*ti5 + ti12*ti4;
 847         cr4 = ti12*tr5 - ti11*tr4;
 848         ci4 = ti12*ti5 - ti11*ti4;
 849         dr3 = cr3 - ci4;
 850         dr4 = cr3 + ci4;
 851         di3 = ci3 + cr4;
 852         di4 = ci3 - cr4;
 853         dr5 = cr2 + ci5;
 854         dr2 = cr2 - ci5;
 855         di5 = ci2 - cr5;
 856         di2 = ci2 + cr5;
 857         ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
 858         ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
 859         ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
 860         ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
 861         ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
 862         ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
 863         ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
 864         ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
 865       }
 866     }
 867   } /* radb5 */
 868
 869
 870 static void radfg(int ido, int ip, int l1, int idl1,
 871       Treal cc[], Treal ch[], const Treal wa[])
 872   {
 873     static const Treal twopi = 6.28318530717959;
 874     int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
 875     Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
 876     arg = twopi / ip;
 877     dcp = cos(arg);
 878     dsp = sin(arg);
 879     ipph = (ip + 1) / 2;
 880     nbd = (ido - 1) / 2;
 881     if (ido != 1) {
 882       for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik];
 883       for (j=1; j<ip; j++)
 884         for (k=0; k<l1; k++)
 885           ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido];
 886       if (nbd <= l1) {
 887         is = -ido;
 888         for (j=1; j<ip; j++) {
 889           is += ido;
 890           idij = is-1;
 891           for (i=2; i<ido; i+=2) {
 892             idij += 2;
 893             for (k=0; k<l1; k++) {
 894               ch[i - 1 + (k + j*l1)*ido] =
 895                   wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
 896               ch[i + (k + j*l1)*ido] =
 897                   wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
 898             }
 899           }
 900         }
 901       } else {
 902         is = -ido;
 903         for (j=1; j<ip; j++) {
 904           is += ido;
 905           for (k=0; k<l1; k++) {
 906             idij = is-1;
 907             for (i=2; i<ido; i+=2) {
 908               idij += 2;
 909               ch[i - 1 + (k + j*l1)*ido] =
 910                   wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido];
 911               ch[i + (k + j*l1)*ido] =
 912                   wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido];
 913             }
 914           }
 915         }
 916       }
 917       if (nbd >= l1) {
 918         for (j=1; j<ipph; j++) {
 919           jc = ip - j;
 920           for (k=0; k<l1; k++) {
 921             for (i=2; i<ido; i+=2) {
 922               cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
 923               cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
 924               cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
 925               cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
 926             }
 927           }
 928         }
 929       } else {
 930         for (j=1; j<ipph; j++) {
 931           jc = ip - j;
 932           for (i=2; i<ido; i+=2) {
 933             for (k=0; k<l1; k++) {
 934               cc[i - 1 + (k + j*l1)*ido] =
 935                   ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
 936               cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido];
 937               cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
 938               cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido];
 939             }
 940           }
 941         }
 942       }
 943     } else {  /* now ido == 1 */
 944       for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
 945     }
 946     for (j=1; j<ipph; j++) {
 947       jc = ip - j;
 948       for (k=0; k<l1; k++) {
 949         cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido];
 950         cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido];
 951       }
 952     }
 953
 954     ar1 = 1;
 955     ai1 = 0;
 956     for (l=1; l<ipph; l++) {
 957       lc = ip - l;
 958       ar1h = dcp*ar1 - dsp*ai1;
 959       ai1 = dcp*ai1 + dsp*ar1;
 960       ar1 = ar1h;
 961       for (ik=0; ik<idl1; ik++) {
 962         ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1];
 963         ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1];
 964       }
 965       dc2 = ar1;
 966       ds2 = ai1;
 967       ar2 = ar1;
 968       ai2 = ai1;
 969       for (j=2; j<ipph; j++) {
 970         jc = ip - j;
 971         ar2h = dc2*ar2 - ds2*ai2;
 972         ai2 = dc2*ai2 + ds2*ar2;
 973         ar2 = ar2h;
 974         for (ik=0; ik<idl1; ik++) {
 975           ch[ik + l*idl1] += ar2*cc[ik + j*idl1];
 976           ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1];
 977         }
 978       }
 979     }
 980     for (j=1; j<ipph; j++)
 981       for (ik=0; ik<idl1; ik++)
 982         ch[ik] += cc[ik + j*idl1];
 983
 984     if (ido >= l1) {
 985       for (k=0; k<l1; k++) {
 986         for (i=0; i<ido; i++) {
 987           ref(cc,i + k*ip*ido) = ch[i + k*ido];
 988         }
 989       }
 990     } else {
 991       for (i=0; i<ido; i++) {
 992         for (k=0; k<l1; k++) {
 993           ref(cc,i + k*ip*ido) = ch[i + k*ido];
 994         }
 995       }
 996     }
 997     for (j=1; j<ipph; j++) {
 998       jc = ip - j;
 999       j2 = 2*j;
1000       for (k=0; k<l1; k++) {
1001         ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
1002             ch[(k + j*l1)*ido];
1003         ref(cc,(j2 + k*ip)*ido) =
1004             ch[(k + jc*l1)*ido];
1005       }
1006     }
1007     if (ido == 1) return;
1008     if (nbd >= l1) {
1009       for (j=1; j<ipph; j++) {
1010         jc = ip - j;
1011         j2 = 2*j;
1012         for (k=0; k<l1; k++) {
1013           for (i=2; i<ido; i+=2) {
1014             ic = ido - i;
1015             ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1016             ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
1017             ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1018             ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
1019           }
1020         }
1021       }
1022     } else {
1023       for (j=1; j<ipph; j++) {
1024         jc = ip - j;
1025         j2 = 2*j;
1026         for (i=2; i<ido; i+=2) {
1027           ic = ido - i;
1028           for (k=0; k<l1; k++) {
1029             ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido];
1030             ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido];
1031             ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido];
1032             ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido];
1033           }
1034         }
1035       }
1036     }
1037   } /* radfg */
1038
1039
1040 static void radbg(int ido, int ip, int l1, int idl1,
1041       Treal cc[], Treal ch[], const Treal wa[])
1042   {
1043     static const Treal twopi = 6.28318530717959;
1044     int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
1045     Treal dc2, ai1, ai2, ar1, ar2, ds2;
1046     int nbd;
1047     Treal dcp, arg, dsp, ar1h, ar2h;
1048     arg = twopi / ip;
1049     dcp = cos(arg);
1050     dsp = sin(arg);
1051     nbd = (ido - 1) / 2;
1052     ipph = (ip + 1) / 2;
1053     if (ido >= l1) {
1054       for (k=0; k<l1; k++) {
1055         for (i=0; i<ido; i++) {
1056           ch[i + k*ido] = ref(cc,i + k*ip*ido);
1057         }
1058       }
1059     } else {
1060       for (i=0; i<ido; i++) {
1061         for (k=0; k<l1; k++) {
1062           ch[i + k*ido] = ref(cc,i + k*ip*ido);
1063         }
1064       }
1065     }
1066     for (j=1; j<ipph; j++) {
1067       jc = ip - j;
1068       j2 = 2*j;
1069       for (k=0; k<l1; k++) {
1070         ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*
1071             ido);
1072         ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
1073       }
1074     }
1075
1076     if (ido != 1) {
1077       if (nbd >= l1) {
1078         for (j=1; j<ipph; j++) {
1079           jc = ip - j;
1080           for (k=0; k<l1; k++) {
1081             for (i=2; i<ido; i+=2) {
1082               ic = ido - i;
1083               ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
1084                   ic - 1 + (2*j - 1 + k*ip)*ido);
1085               ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
1086                   ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
1087               ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
1088                   + (2*j - 1 + k*ip)*ido);
1089               ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
1090                   + (2*j - 1 + k*ip)*ido);
1091             }
1092           }
1093         }
1094       } else {
1095         for (j=1; j<ipph; j++) {
1096           jc = ip - j;
1097           for (i=2; i<ido; i+=2) {
1098             ic = ido - i;
1099             for (k=0; k<l1; k++) {
1100               ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,
1101                   ic - 1 + (2*j - 1 + k*ip)*ido);
1102               ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) -
1103                   ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
1104               ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic
1105                   + (2*j - 1 + k*ip)*ido);
1106               ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic
1107                   + (2*j - 1 + k*ip)*ido);
1108             }
1109           }
1110         }
1111       }
1112     }
1113
1114     ar1 = 1;
1115     ai1 = 0;
1116     for (l=1; l<ipph; l++) {
1117       lc = ip - l;
1118       ar1h = dcp*ar1 - dsp*ai1;
1119       ai1 = dcp*ai1 + dsp*ar1;
1120       ar1 = ar1h;
1121       for (ik=0; ik<idl1; ik++) {
1122         cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1];
1123         cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1];
1124       }
1125       dc2 = ar1;
1126       ds2 = ai1;
1127       ar2 = ar1;
1128       ai2 = ai1;
1129       for (j=2; j<ipph; j++) {
1130         jc = ip - j;
1131         ar2h = dc2*ar2 - ds2*ai2;
1132         ai2 = dc2*ai2 + ds2*ar2;
1133         ar2 = ar2h;
1134         for (ik=0; ik<idl1; ik++) {
1135           cc[ik + l*idl1] += ar2*ch[ik + j*idl1];
1136           cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1];
1137         }
1138       }
1139     }
1140     for (j=1; j<ipph; j++) {
1141       for (ik=0; ik<idl1; ik++) {
1142         ch[ik] += ch[ik + j*idl1];
1143       }
1144     }
1145     for (j=1; j<ipph; j++) {
1146       jc = ip - j;
1147       for (k=0; k<l1; k++) {
1148         ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido];
1149         ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido];
1150       }
1151     }
1152
1153     if (ido == 1) return;
1154     if (nbd >= l1) {
1155       for (j=1; j<ipph; j++) {
1156         jc = ip - j;
1157         for (k=0; k<l1; k++) {
1158           for (i=2; i<ido; i+=2) {
1159             ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
1160             ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido];
1161             ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
1162             ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
1163           }
1164         }
1165       }
1166     } else {
1167       for (j=1; j<ipph; j++) {
1168         jc = ip - j;
1169         for (i=2; i<ido; i+=2) {
1170           for (k=0; k<l1; k++) {
1171             ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido];
1172             ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido];
1173             ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido];
1174             ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido];
1175           }
1176         }
1177       }
1178     }
1179     for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik];
1180     for (j=1; j<ip; j++)
1181       for (k=0; k<l1; k++)
1182         cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido];
1183     if (nbd <= l1) {
1184       is = -ido;
1185       for (j=1; j<ip; j++) {
1186         is += ido;
1187         idij = is-1;
1188         for (i=2; i<ido; i+=2) {
1189           idij += 2;
1190           for (k=0; k<l1; k++) {
1191             cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
1192                 ch[i + (k + j*l1)*ido];
1193             cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
1194           }
1195         }
1196       }
1197     } else {
1198       is = -ido;
1199       for (j=1; j<ip; j++) {
1200         is += ido;
1201         for (k=0; k<l1; k++) {
1202           idij = is - 1;
1203           for (i=2; i<ido; i+=2) {
1204             idij += 2;
1205             cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]*
1206                 ch[i + (k + j*l1)*ido];
1207             cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido];
1208           }
1209         }
1210       }
1211     }
1212   } /* radbg */
1213
1214   /* ----------------------------------------------------------------------
1215 cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
1216 ---------------------------------------------------------------------- */
1217
1218 void fftpack_cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign)
1219   {
1220     int idot, i;
1221     int k1, l1, l2;
1222     int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
1223     Treal *cinput, *coutput;
1224     nf = ifac[1];
1225     na = 0;
1226     l1 = 1;
1227     iw = 0;
1228     for (k1=2; k1<=nf+1; k1++) {
1229       ip = ifac[k1];
1230       l2 = ip*l1;
1231       ido = n / l2;
1232       idot = ido + ido;
1233       idl1 = idot*l1;
1234       if (na) {
1235         cinput = ch;
1236         coutput = c;
1237       } else {
1238         cinput = c;
1239         coutput = ch;
1240       }
1241       switch (ip) {
1242       case 4:
1243         ix2 = iw + idot;
1244         ix3 = ix2 + idot;
1245         passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign);
1246         na = !na;
1247         break;
1248       case 2:
1249         passf2(idot, l1, cinput, coutput, &wa[iw], isign);
1250         na = !na;
1251         break;
1252       case 3:
1253         ix2 = iw + idot;
1254         passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign);
1255         na = !na;
1256         break;
1257       case 5:
1258         ix2 = iw + idot;
1259         ix3 = ix2 + idot;
1260         ix4 = ix3 + idot;
1261         passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
1262         na = !na;
1263         break;
1264       default:
1265         passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign);
1266         if (nac != 0) na = !na;
1267       }
1268       l1 = l2;
1269       iw += (ip - 1)*idot;
1270     }
1271     if (na == 0) return;
1272     for (i=0; i<2*n; i++) c[i] = ch[i];
1273   } /* cfftf1 */
1274
1275
1276 void fftpack_cfftf(int n, Treal c[], Treal wsave[])
1277   {
1278     int iw1, iw2;
1279     if (n == 1) return;
1280     iw1 = 2*n;
1281     iw2 = iw1 + 2*n;
1282     fftpack_cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1);
1283   } /* cfftf */
1284
1285
1286 void fftpack_cfftb(int n, Treal c[], Treal wsave[])
1287   {
1288     int iw1, iw2;
1289     if (n == 1) return;
1290     iw1 = 2*n;
1291     iw2 = iw1 + 2*n;
1292     fftpack_cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1);
1293   } /* cfftb */
1294
1295
1296 static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL])
1297   /* Factorize n in factors in ntryh and rest. On exit,
1298 ifac[0] contains n and ifac[1] contains number of factors,
1299 the factors start from ifac[2]. */
1300   {
1301     int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
1302 startloop:
1303     if (j < NSPECIAL)
1304       ntry = ntryh[j];
1305     else
1306       ntry+= 2;
1307     j++;
1308     do {
1309       nq = nl / ntry;
1310       nr = nl - ntry*nq;
1311       if (nr != 0) goto startloop;
1312       nf++;
1313       ifac[nf + 1] = ntry;
1314       nl = nq;
1315       if (ntry == 2 && nf != 1) {
1316         for (i=2; i<=nf; i++) {
1317           ib = nf - i + 2;
1318           ifac[ib + 1] = ifac[ib];
1319         }
1320         ifac[2] = 2;
1321       }
1322     } while (nl != 1);
1323     ifac[0] = n;
1324     ifac[1] = nf;
1325   }
1326
1327
1328 void fftpack_cffti1(int n, Treal wa[], int ifac[MAXFAC+2])
1329   {
1330     static const Treal twopi = 6.28318530717959;
1331     Treal arg, argh, argld, fi;
1332     int idot, i, j;
1333     int i1, k1, l1, l2;
1334     int ld, ii, nf, ip;
1335     int ido, ipm;
1336
1337     static const int ntryh[NSPECIAL] = {
1338       3,4,2,5    }; /* Do not change the order of these. */
1339
1340     factorize(n,ifac,ntryh);
1341     nf = ifac[1];
1342     argh = twopi/(Treal)n;
1343     i = 1;
1344     l1 = 1;
1345     for (k1=1; k1<=nf; k1++) {
1346       ip = ifac[k1+1];
1347       ld = 0;
1348       l2 = l1*ip;
1349       ido = n / l2;
1350       idot = ido + ido + 2;
1351       ipm = ip - 1;
1352       for (j=1; j<=ipm; j++) {
1353         i1 = i;
1354         wa[i-1] = 1;
1355         wa[i] = 0;
1356         ld += l1;
1357         fi = 0;
1358         argld = ld*argh;
1359         for (ii=4; ii<=idot; ii+=2) {
1360           i+= 2;
1361           fi+= 1;
1362           arg = fi*argld;
1363           wa[i-1] = cos(arg);
1364           wa[i] = sin(arg);
1365         }
1366         if (ip > 5) {
1367           wa[i1-1] = wa[i-1];
1368           wa[i1] = wa[i];
1369         }
1370       }
1371       l1 = l2;
1372     }
1373   } /* cffti1 */
1374
1375
1376 void fftpack_cffti(int n, Treal wsave[])
1377  {
1378     int iw1, iw2;
1379     if (n == 1) return;
1380     iw1 = 2*n;
1381     iw2 = iw1 + 2*n;
1382     fftpack_cffti1(n, wsave+iw1, (int*)(wsave+iw2));
1383   } /* cffti */
1384
1385   /* ----------------------------------------------------------------------
1386 rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Treal FFTs.
1387 ---------------------------------------------------------------------- */
1388
1389 void fftpack_rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
1390   {
1391     int i;
1392     int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
1393     Treal *cinput, *coutput;
1394     nf = ifac[1];
1395     na = 1;
1396     l2 = n;
1397     iw = n-1;
1398     for (k1 = 1; k1 <= nf; ++k1) {
1399       kh = nf - k1;
1400       ip = ifac[kh + 2];
1401       l1 = l2 / ip;
1402       ido = n / l2;
1403       idl1 = ido*l1;
1404       iw -= (ip - 1)*ido;
1405       na = !na;
1406       if (na) {
1407         cinput = ch;
1408         coutput = c;
1409       } else {
1410         cinput = c;
1411         coutput = ch;
1412       }
1413       switch (ip) {
1414       case 4:
1415         ix2 = iw + ido;
1416         ix3 = ix2 + ido;
1417         radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
1418         break;
1419       case 2:
1420         radf2(ido, l1, cinput, coutput, &wa[iw]);
1421         break;
1422       case 3:
1423         ix2 = iw + ido;
1424         radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
1425         break;
1426       case 5:
1427         ix2 = iw + ido;
1428         ix3 = ix2 + ido;
1429         ix4 = ix3 + ido;
1430         radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
1431         break;
1432       default:
1433         if (ido == 1)
1434           na = !na;
1435         if (na == 0) {
1436           radfg(ido, ip, l1, idl1, c, ch, &wa[iw]);
1437           na = 1;
1438         } else {
1439           radfg(ido, ip, l1, idl1, ch, c, &wa[iw]);
1440           na = 0;
1441         }
1442       }
1443       l2 = l1;
1444     }
1445     if (na == 1) return;
1446     for (i = 0; i < n; i++) c[i] = ch[i];
1447   } /* rfftf1 */
1448
1449
1450 void fftpack_rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2])
1451   {
1452     int i;
1453     int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
1454     Treal *cinput, *coutput;
1455     nf = ifac[1];
1456     na = 0;
1457     l1 = 1;
1458     iw = 0;
1459     for (k1=1; k1<=nf; k1++) {
1460       ip = ifac[k1 + 1];
1461       l2 = ip*l1;
1462       ido = n / l2;
1463       idl1 = ido*l1;
1464       if (na) {
1465         cinput = ch;
1466         coutput = c;
1467       } else {
1468         cinput = c;
1469         coutput = ch;
1470       }
1471       switch (ip) {
1472       case 4:
1473         ix2 = iw + ido;
1474         ix3 = ix2 + ido;
1475         radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]);
1476         na = !na;
1477         break;
1478       case 2:
1479         radb2(ido, l1, cinput, coutput, &wa[iw]);
1480         na = !na;
1481         break;
1482       case 3:
1483         ix2 = iw + ido;
1484         radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]);
1485         na = !na;
1486         break;
1487       case 5:
1488         ix2 = iw + ido;
1489         ix3 = ix2 + ido;
1490         ix4 = ix3 + ido;
1491         radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]);
1492         na = !na;
1493         break;
1494       default:
1495         radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]);
1496         if (ido == 1) na = !na;
1497       }
1498       l1 = l2;
1499       iw += (ip - 1)*ido;
1500     }
1501     if (na == 0) return;
1502     for (i=0; i<n; i++) c[i] = ch[i];
1503   } /* rfftb1 */
1504
1505
1506 void fftpack_rfftf(int n, Treal r[], Treal wsave[])
1507   {
1508     if (n == 1) return;
1509     fftpack_rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
1510   } /* rfftf */
1511
1512
1513 void fftpack_rfftb(int n, Treal r[], Treal wsave[])
1514   {
1515     if (n == 1) return;
1516     fftpack_rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n));
1517   } /* rfftb */
1518
1519
1520 void fftpack_rffti1(int n, Treal wa[], int ifac[MAXFAC+2])
1521   {
1522     static const Treal twopi = 6.28318530717959;
1523     Treal arg, argh, argld, fi;
1524     int i, j;
1525     int k1, l1, l2;
1526     int ld, ii, nf, ip, is;
1527     int ido, ipm, nfm1;
1528     static const int ntryh[NSPECIAL] = {
1529       4,2,3,5    }; /* Do not change the order of these. */
1530     factorize(n,ifac,ntryh);
1531     nf = ifac[1];
1532     argh = twopi / n;
1533     is = 0;
1534     nfm1 = nf - 1;
1535     l1 = 1;
1536     if (nfm1 == 0) return;
1537     for (k1 = 1; k1 <= nfm1; k1++) {
1538       ip = ifac[k1 + 1];
1539       ld = 0;
1540       l2 = l1*ip;
1541       ido = n / l2;
1542       ipm = ip - 1;
1543       for (j = 1; j <= ipm; ++j) {
1544         ld += l1;
1545         i = is;
1546         argld = (Treal) ld*argh;
1547         fi = 0;
1548         for (ii = 3; ii <= ido; ii += 2) {
1549           i += 2;
1550           fi += 1;
1551           arg = fi*argld;
1552           wa[i - 2] = cos(arg);
1553           wa[i - 1] = sin(arg);
1554         }
1555         is += ido;
1556       }
1557       l1 = l2;
1558     }
1559   } /* rffti1 */
1560
1561
1562 void fftpack_rffti(int n, Treal wsave[])
1563   {
1564     if (n == 1) return;
1565     fftpack_rffti1(n, wsave+n, (int*)(wsave+2*n));
1566   } /* rffti */
1567
1568 #ifdef __cplusplus
1569 }
1570 #endif