src/external/fftpack/fftpack.c

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