src/external/fftpack/fftpack.cpp

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