lib/cfft.c

   1 /* from fftpkg package in cmlib and netlib -- translated by f2c and modified */
   2
   3 #include "xmath.h"
   4
   5 /*
   6   Internal Routines
   7 */
   8
   9 static int cffti1_(int *n, double *wa, double *ifac)
  10 {
  11   /* Initialized data */
  12   static int ntryh[4] = { 3,4,2,5 };
  13
  14   /* System generated locals */
  15   long i_1, i_2, i_3;
  16
  17   /* Local variables */
  18   double argh;
  19   long ntry;
  20   int i, j;
  21   long idot;
  22   double argld;
  23   int i1, ib;
  24   long k1, l1, l2;
  25   double fi;
  26   int ii, nf;
  27   long ip, ld, nl, nq, nr;
  28   double arg;
  29   long ido, ipm;
  30   double tpi;
  31
  32   /* Parameter adjustments */
  33   --wa;
  34   --ifac;
  35
  36   /* Function Body */
  37   nl = *n;
  38   nf = 0;
  39   j = 0;
  40
  41  L101:
  42   ++j;
  43   if (j - 4 <= 0) ntry = ntryh[j - 1];
  44   else ntry += 2;
  45  L104:
  46   nq = nl / ntry;
  47   nr = nl - ntry * nq;
  48   if (nr != 0) goto L101;
  49   ++nf;
  50   ifac[nf + 2] = ntry;
  51   nl = nq;
  52   if (ntry == 2 && nf != 1) {
  53     i_1 = nf;
  54     for (i = 2; i <= i_1; ++i) {
  55       ib = nf - i + 2;
  56       ifac[ib + 2] = ifac[ib + 1];
  57     }
  58     ifac[3] = 2;
  59   }
  60   if (nl != 1) goto L104;
  61
  62   ifac[1] = *n;
  63   ifac[2] = nf;
  64   tpi = 6.28318530717959;
  65   argh = tpi / (double) (*n);
  66   i = 2;
  67   l1 = 1;
  68   i_1 = nf;
  69   for (k1 = 1; k1 <= i_1; ++k1) {
  70     ip = ifac[k1 + 2];
  71     ld = 0;
  72     l2 = l1 * ip;
  73     ido = *n / l2;
  74     idot = ido + ido + 2;
  75     ipm = ip - 1;
  76     i_2 = ipm;
  77     for (j = 1; j <= i_2; ++j) {
  78       i1 = i;
  79       wa[i - 1] = 1.0;
  80       wa[i] = 0.0;
  81       ld += l1;
  82       fi = 0.0;
  83       argld = (double) ld * argh;
  84       i_3 = idot;
  85       for (ii = 4; ii <= i_3; ii += 2) {
  86         i += 2;
  87         fi += 1.0;
  88         arg = fi * argld;
  89         wa[i - 1] = cos(arg);
  90         wa[i] = sin(arg);
  91       }
  92       if (ip > 5) {
  93         wa[i1 - 1] = wa[i - 1];
  94         wa[i1] = wa[i];
  95       }
  96     }
  97     l1 = l2;
  98   }
  99   return 0;
 100 }
 101
 102 static int pass_(int *nac, long *ido, long *ip, long *l1, long *idl1,
 103                  double *cc, double *c1, double *c2, double *ch,
 104                  double *ch2, double *wa, int *isign)
 105 {
 106   /* System generated locals */
 107   long cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_offset,
 108        ch2_offset, c1_dim1, c1_dim2, c1_offset, ch2_dim1, c2_offset,
 109        c2_dim1, i_1, i_2, i_3;
 110
 111   /* Local variables */
 112   int i, j, k, l,  ik;
 113   long idij, idlj, idot, ipph,   lc, jc,    idj, idl, inc, idp;
 114   double wai, war;
 115   long ipp2;
 116
 117   /* Parameter adjustments */
 118   cc_dim1 = *ido;
 119   cc_dim2 = *ip;
 120   cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
 121   cc -= cc_offset;
 122   c1_dim1 = *ido;
 123   c1_dim2 = *l1;
 124   c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
 125   c1 -= c1_offset;
 126   c2_dim1 = *idl1;
 127   c2_offset = c2_dim1 + 1;
 128   c2 -= c2_offset;
 129   ch_dim1 = *ido;
 130   ch_dim2 = *l1;
 131   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 132   ch -= ch_offset;
 133   ch2_dim1 = *idl1;
 134   ch2_offset = ch2_dim1 + 1;
 135   ch2 -= ch2_offset;
 136   --wa;
 137
 138   /* Function Body */
 139   idot = *ido / 2;
 140   ipp2 = *ip + 2;
 141   ipph = (*ip + 1) / 2;
 142   idp = *ip * *ido;
 143
 144   if (*ido >= *l1) {
 145     i_1 = ipph;
 146     for (j = 2; j <= i_1; ++j) {
 147       jc = ipp2 - j;
 148       i_2 = *l1;
 149       for (k = 1; k <= i_2; ++k) {
 150         i_3 = *ido;
 151         for (i = 1; i <= i_3; ++i) {
 152           ch[i + (k + j * ch_dim2) * ch_dim1] =
 153             cc[i + (j + k * cc_dim2) * cc_dim1]
 154               + cc[i + (jc + k * cc_dim2) * cc_dim1];
 155           ch[i + (k + jc * ch_dim2) * ch_dim1] =
 156             cc[i + (j + k * cc_dim2) * cc_dim1]
 157               - cc[i + (jc + k * cc_dim2) * cc_dim1];
 158         }
 159       }
 160     }
 161     i_1 = *l1;
 162     for (k = 1; k <= i_1; ++k) {
 163       i_2 = *ido;
 164       for (i = 1; i <= i_2; ++i) {
 165         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
 166       }
 167     }
 168   }
 169   else {
 170     i_1 = ipph;
 171     for (j = 2; j <= i_1; ++j) {
 172       jc = ipp2 - j;
 173       i_2 = *ido;
 174       for (i = 1; i <= i_2; ++i) {
 175         i_3 = *l1;
 176         for (k = 1; k <= i_3; ++k) {
 177           ch[i + (k + j * ch_dim2) * ch_dim1] =
 178             cc[i + (j + k * cc_dim2) * cc_dim1]
 179               + cc[i + (jc + k * cc_dim2) * cc_dim1];
 180           ch[i + (k + jc * ch_dim2) * ch_dim1] =
 181             cc[i + (j + k * cc_dim2) * cc_dim1]
 182               - cc[i + (jc + k * cc_dim2) * cc_dim1];
 183         }
 184       }
 185     }
 186     i_1 = *ido;
 187     for (i = 1; i <= i_1; ++i) {
 188       i_2 = *l1;
 189       for (k = 1; k <= i_2; ++k) {
 190         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
 191       }
 192     }
 193   }
 194   idl = 2 - *ido;
 195   inc = 0;
 196   i_1 = ipph;
 197   for (l = 2; l <= i_1; ++l) {
 198     lc = ipp2 - l;
 199     idl += *ido;
 200     i_2 = *idl1;
 201     for (ik = 1; ik <= i_2; ++ik) {
 202       c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1]
 203         * ch2[ik + (ch2_dim1 << 1)];
 204       c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1];
 205     }
 206     idlj = idl;
 207     inc += *ido;
 208     i_2 = ipph;
 209     for (j = 3; j <= i_2; ++j) {
 210       jc = ipp2 - j;
 211       idlj += inc;
 212       if (idlj > idp) {
 213         idlj -= idp;
 214       }
 215       war = wa[idlj - 1];
 216       wai = wa[idlj];
 217       i_3 = *idl1;
 218       for (ik = 1; ik <= i_3; ++ik) {
 219         c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
 220         c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1];
 221       }
 222     }
 223   }
 224   i_1 = ipph;
 225   for (j = 2; j <= i_1; ++j) {
 226     i_2 = *idl1;
 227     for (ik = 1; ik <= i_2; ++ik) {
 228       ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
 229     }
 230   }
 231   i_1 = ipph;
 232   for (j = 2; j <= i_1; ++j) {
 233     jc = ipp2 - j;
 234     i_2 = *idl1;
 235     for (ik = 2; ik <= i_2; ik += 2) {
 236       ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
 237         - c2[ik + jc * c2_dim1];
 238       ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
 239         + c2[ik + jc * c2_dim1];
 240       ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1]
 241         + c2[ik - 1 + jc * c2_dim1];
 242       ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1]
 243         - c2[ik - 1 + jc * c2_dim1];
 244     }
 245   }
 246   *nac = 1;
 247   if (*ido != 2) {
 248     *nac = 0;
 249     i_1 = *idl1;
 250     for (ik = 1; ik <= i_1; ++ik) {
 251       c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
 252     }
 253     i_1 = *ip;
 254     for (j = 2; j <= i_1; ++j) {
 255       i_2 = *l1;
 256       for (k = 1; k <= i_2; ++k) {
 257         c1[(k + j * c1_dim2) * c1_dim1 + 1] =
 258           ch[(k + j * ch_dim2) * ch_dim1 + 1];
 259         c1[(k + j * c1_dim2) * c1_dim1 + 2] =
 260           ch[(k + j * ch_dim2) * ch_dim1 + 2];
 261       }
 262     }
 263     if (idot <= *l1) {
 264       idij = 0;
 265       i_1 = *ip;
 266       for (j = 2; j <= i_1; ++j) {
 267         idij += 2;
 268         i_2 = *ido;
 269         for (i = 4; i <= i_2; i += 2) {
 270           idij += 2;
 271           i_3 = *l1;
 272           for (k = 1; k <= i_3; ++k) {
 273             c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 274               * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij]
 275                 * ch[i + (k + j * ch_dim2) * ch_dim1];
 276             c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 277               * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
 278                 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
 279           }
 280         }
 281       }
 282     }
 283     else {
 284       idj = 2 - *ido;
 285       i_1 = *ip;
 286       for (j = 2; j <= i_1; ++j) {
 287         idj += *ido;
 288         i_2 = *l1;
 289         for (k = 1; k <= i_2; ++k) {
 290           idij = idj;
 291           i_3 = *ido;
 292           for (i = 4; i <= i_3; i += 2) {
 293             idij += 2;
 294             c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 295               * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]
 296                 + *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1];
 297             c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 298               * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
 299                 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
 300           }
 301         }
 302       }
 303     }
 304   }
 305   return 0;
 306 }
 307
 308 static int pass2_(ido, l1, cc, ch, wa1, isign)
 309      int *ido, *l1;
 310      double *cc, *ch, *wa1;
 311      int *isign;
 312 {
 313   /* System generated locals */
 314   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 315
 316   /* Local variables */
 317   int i, k;
 318   double ti2, tr2;
 319
 320   /* Parameter adjustments */
 321   cc_dim1 = *ido;
 322   cc_offset = cc_dim1 * 3 + 1;
 323   cc -= cc_offset;
 324   ch_dim1 = *ido;
 325   ch_dim2 = *l1;
 326   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 327   ch -= ch_offset;
 328   --wa1;
 329
 330   /* Function Body */
 331   if (*ido <= 2) {
 332     i_1 = *l1;
 333     for (k = 1; k <= i_1; ++k) {
 334       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
 335         + cc[((k << 1) + 2) * cc_dim1 + 1];
 336       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
 337         - cc[((k << 1) + 2) * cc_dim1 + 1];
 338       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
 339         + cc[((k << 1) + 2) * cc_dim1 + 2];
 340       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
 341         - cc[((k << 1) + 2) * cc_dim1 + 2];
 342     }
 343   }
 344   else {
 345     i_1 = *l1;
 346     for (k = 1; k <= i_1; ++k) {
 347       i_2 = *ido;
 348       for (i = 2; i <= i_2; i += 2) {
 349         ch[i - 1 + (k + ch_dim2) * ch_dim1]
 350           = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
 351             + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
 352         tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
 353           - cc[i - 1 + ((k << 1) + 2) * cc_dim1];
 354         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
 355           + cc[i + ((k << 1) + 2) * cc_dim1];
 356         ti2 = cc[i + ((k << 1) + 1) * cc_dim1]
 357           - cc[i + ((k << 1) + 2) * cc_dim1];
 358         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2
 359           - *isign * wa1[i] * tr2;
 360         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2
 361           + *isign * wa1[i] * ti2;
 362       }
 363     }
 364   }
 365   return 0;
 366 }
 367
 368 static int pass3_(ido, l1, cc, ch, wa1, wa2, isign)
 369      int *ido, *l1;
 370      double *cc, *ch, *wa1, *wa2;
 371      int *isign;
 372 {
 373   /* System generated locals */
 374   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 375
 376   /* Local variables */
 377   double taui, taur;
 378   int i, k;
 379   double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
 380
 381   /* Parameter adjustments */
 382   cc_dim1 = *ido;
 383   cc_offset = (cc_dim1 << 2) + 1;
 384   cc -= cc_offset;
 385   ch_dim1 = *ido;
 386   ch_dim2 = *l1;
 387   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 388   ch -= ch_offset;
 389   --wa1;
 390   --wa2;
 391
 392   /* Function Body */
 393   taur = -.5;
 394   taui = -(*isign) * .866025403784439;
 395   if (*ido == 2) {
 396     i_1 = *l1;
 397     for (k = 1; k <= i_1; ++k) {
 398       tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
 399       cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
 400       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
 401       ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
 402       ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
 403       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
 404       cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1]
 405                     - cc[(k * 3 + 3) * cc_dim1 + 1]);
 406       ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2]
 407                     - cc[(k * 3 + 3) * cc_dim1 + 2]);
 408       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
 409       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
 410       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
 411       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
 412     }
 413   }
 414   else {
 415     i_1 = *l1;
 416     for (k = 1; k <= i_1; ++k) {
 417       i_2 = *ido;
 418       for (i = 2; i <= i_2; i += 2) {
 419         tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1]
 420           + cc[i - 1 + (k * 3 + 3) * cc_dim1];
 421         cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
 422         ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1)
 423                                                  * cc_dim1] + tr2;
 424         ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1];
 425         ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
 426         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2;
 427         cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1]
 428                       - cc[i - 1 + (k * 3 + 3) * cc_dim1]);
 429         ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1]
 430                       - cc[i + (k * 3 + 3) * cc_dim1]);
 431         dr2 = cr2 - ci3;
 432         dr3 = cr2 + ci3;
 433         di2 = ci2 + cr3;
 434         di3 = ci2 - cr3;
 435         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
 436           - *isign * wa1[i] * dr2;
 437         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
 438           + *isign * wa1[i] * di2;
 439         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
 440           - *isign * wa2[i] * dr3;
 441         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
 442           + *isign * wa2[i] * di3;
 443       }
 444     }
 445   }
 446   return 0;
 447 }
 448
 449 static int pass4_(ido, l1, cc, ch, wa1, wa2, wa3, isign)
 450      int *ido, *l1;
 451      double *cc, *ch, *wa1, *wa2, *wa3;
 452      int *isign;
 453 {
 454   /* System generated locals */
 455   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 456
 457   /* Local variables */
 458   int i, k;
 459   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
 460
 461   /* Parameter adjustments */
 462   cc_dim1 = *ido;
 463   cc_offset = cc_dim1 * 5 + 1;
 464   cc -= cc_offset;
 465   ch_dim1 = *ido;
 466   ch_dim2 = *l1;
 467   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 468   ch -= ch_offset;
 469   --wa1;
 470   --wa2;
 471   --wa3;
 472
 473   /* Function Body */
 474   if (*ido == 2) {
 475     i_1 = *l1;
 476     for (k = 1; k <= i_1; ++k) {
 477       ti1 = cc[((k << 2) + 1) * cc_dim1 + 2]
 478         - cc[((k << 2) + 3) * cc_dim1 + 2];
 479       ti2 = cc[((k << 2) + 1) * cc_dim1 + 2]
 480         + cc[((k << 2) + 3) * cc_dim1 + 2];
 481       tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2]
 482                       - cc[((k << 2) + 4) * cc_dim1 + 2]);
 483       ti3 = cc[((k << 2) + 2) * cc_dim1 + 2]
 484         + cc[((k << 2) + 4) * cc_dim1 + 2];
 485       tr1 = cc[((k << 2) + 1) * cc_dim1 + 1]
 486         - cc[((k << 2) + 3) * cc_dim1 + 1];
 487       tr2 = cc[((k << 2) + 1) * cc_dim1 + 1]
 488         + cc[((k << 2) + 3) * cc_dim1 + 1];
 489       ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1]
 490                       - cc[((k << 2) + 2) * cc_dim1 + 1]);
 491       tr3 = cc[((k << 2) + 2) * cc_dim1 + 1]
 492         + cc[((k << 2) + 4) * cc_dim1 + 1];
 493       ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
 494       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
 495       ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
 496       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
 497       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
 498       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
 499       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
 500       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
 501     }
 502   }
 503   else {
 504     i_1 = *l1;
 505     for (k = 1; k <= i_1; ++k) {
 506       i_2 = *ido;
 507       for (i = 2; i <= i_2; i += 2) {
 508         ti1 = cc[i + ((k << 2) + 1) * cc_dim1]
 509           - cc[i + ((k << 2) + 3) * cc_dim1];
 510         ti2 = cc[i + ((k << 2) + 1) * cc_dim1]
 511           + cc[i + ((k << 2) + 3) * cc_dim1];
 512         ti3 = cc[i + ((k << 2) + 2) * cc_dim1]
 513           + cc[i + ((k << 2) + 4) * cc_dim1];
 514         tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1]
 515                         - cc[i + ((k << 2) + 4) * cc_dim1]);
 516         tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
 517           - cc[i - 1 + ((k << 2) + 3) * cc_dim1];
 518         tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
 519           + cc[i - 1 + ((k << 2) + 3) * cc_dim1];
 520         ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1]
 521                         - cc[i - 1 + ((k << 2) + 2) * cc_dim1]);
 522         tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1]
 523           + cc[i - 1 + ((k << 2) + 4) * cc_dim1];
 524         ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
 525         cr3 = tr2 - tr3;
 526         ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
 527         ci3 = ti2 - ti3;
 528         cr2 = tr1 + tr4;
 529         cr4 = tr1 - tr4;
 530         ci2 = ti1 + ti4;
 531         ci4 = ti1 - ti4;
 532         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2
 533           + *isign * wa1[i] * ci2;
 534         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2
 535           - *isign * wa1[i] * cr2;
 536         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3
 537           + *isign * wa2[i] * ci3;
 538         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3
 539           - *isign * wa2[i] * cr3;
 540         ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4
 541           + *isign * wa3[i] * ci4;
 542         ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4
 543           - *isign * wa3[i] * cr4;
 544       }
 545     }
 546   }
 547   return 0;
 548 }
 549
 550 static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign)
 551      int *ido, *l1;
 552      double *cc, *ch, *wa1, *wa2, *wa3, *wa4;
 553      int *isign;
 554 {
 555   /* System generated locals */
 556   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 557
 558   /* Local variables */
 559   int i, k;
 560   double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
 561          ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11,
 562          tr12;
 563
 564   /* Parameter adjustments */
 565   cc_dim1 = *ido;
 566   cc_offset = cc_dim1 * 6 + 1;
 567   cc -= cc_offset;
 568   ch_dim1 = *ido;
 569   ch_dim2 = *l1;
 570   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 571   ch -= ch_offset;
 572   --wa1;
 573   --wa2;
 574   --wa3;
 575   --wa4;
 576
 577   /* Function Body */
 578   tr11 = .309016994374947;
 579   ti11 = -(*isign) * .951056516295154;
 580   tr12 = -.809016994374947;
 581   ti12 = -(*isign) * .587785252292473;
 582   if (*ido == 2) {
 583     i_1 = *l1;
 584     for (k = 1; k <= i_1; ++k) {
 585       ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
 586       ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
 587       ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
 588       ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
 589       tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
 590       tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
 591       tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
 592       tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
 593       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1]
 594         + tr2 + tr3;
 595       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2]
 596         + ti2 + ti3;
 597       cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
 598       ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
 599       cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
 600       ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
 601       cr5 = ti11 * tr5 + ti12 * tr4;
 602       ci5 = ti11 * ti5 + ti12 * ti4;
 603       cr4 = ti12 * tr5 - ti11 * tr4;
 604       ci4 = ti12 * ti5 - ti11 * ti4;
 605       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
 606       ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
 607       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
 608       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
 609       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
 610       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
 611       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
 612       ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
 613     }
 614   }
 615   else {
 616     i_1 = *l1;
 617     for (k = 1; k <= i_1; ++k) {
 618       i_2 = *ido;
 619       for (i = 2; i <= i_2; i += 2) {
 620         ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1];
 621         ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1];
 622         ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1];
 623         ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1];
 624         tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
 625           - cc[i - 1 + (k * 5 + 5) * cc_dim1];
 626         tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
 627           + cc[i - 1 + (k * 5 + 5) * cc_dim1];
 628         tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
 629           - cc[i - 1 + (k * 5 + 4) * cc_dim1];
 630         tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
 631           + cc[i - 1 + (k * 5 + 4) * cc_dim1];
 632         ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1)
 633                                                  * cc_dim1] + tr2 + tr3;
 634         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1)
 635                                              * cc_dim1] + ti2 + ti3;
 636         cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
 637         ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
 638
 639         cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
 640         ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
 641
 642         cr5 = ti11 * tr5 + ti12 * tr4;
 643         ci5 = ti11 * ti5 + ti12 * ti4;
 644         cr4 = ti12 * tr5 - ti11 * tr4;
 645         ci4 = ti12 * ti5 - ti11 * ti4;
 646         dr3 = cr3 - ci4;
 647         dr4 = cr3 + ci4;
 648         di3 = ci3 + cr4;
 649         di4 = ci3 - cr4;
 650         dr5 = cr2 + ci5;
 651         dr2 = cr2 - ci5;
 652         di5 = ci2 - cr5;
 653         di2 = ci2 + cr5;
 654         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
 655           + *isign * wa1[i] * di2;
 656         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
 657           - *isign * wa1[i] * dr2;
 658         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
 659           + *isign * wa2[i] * di3;
 660         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
 661           - *isign * wa2[i] * dr3;
 662         ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4
 663           + *isign * wa3[i] * di4;
 664         ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4
 665           - *isign * wa3[i] * dr4;
 666         ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5
 667           + *isign * wa4[i] * di5;
 668         ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5
 669           - *isign * wa4[i] * dr5;
 670       }
 671     }
 672   }
 673   return 0;
 674 }
 675
 676
 677 /*static */
 678 int
 679 cfft1_(int *n, double *c, double *ch, double *wa,
 680        double /* int??*/  *ifac, int *isign)
 681 {
 682   /* System generated locals */
 683   long i_1;
 684
 685   /* Local variables */
 686   long idot;
 687   int i, k1, n2, na, nac;
 688   long l1, l2, ido, idl1, iw, ix2, ix3, ix4;
 689   long nf, ip;
 690
 691   /* Parameter adjustments */
 692   --c;
 693   --ch;
 694   --wa;
 695   --ifac;
 696
 697   /* Function Body */
 698   nf = ifac[2];
 699   na = 0;
 700   l1 = 1;
 701   iw = 1;
 702   i_1 = nf;
 703   for (k1 = 1; k1 <= i_1; ++k1) {
 704     ip = ifac[k1 + 2];
 705     l2 = ip * l1;
 706     ido = *n / l2;
 707     idot = ido + ido;
 708     idl1 = idot * l1;
 709     if (ip == 4) {
 710       ix2 = iw + idot;
 711       ix3 = ix2 + idot;
 712       if (na == 0) {
 713         pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
 714       }
 715       else {
 716         pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
 717       }
 718       na = 1 - na;
 719     }
 720     else if (ip == 2) {
 721       if (na == 0) {
 722         pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign);
 723       }
 724       else {
 725         pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign);
 726       }
 727       na = 1 - na;
 728     }
 729     else if (ip == 3) {
 730       ix2 = iw + idot;
 731       if (na == 0) {
 732         pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign);
 733       }
 734       else {
 735         pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign);
 736       }
 737       na = 1 - na;
 738     }
 739     else if (ip == 5) {
 740       ix2 = iw + idot;
 741       ix3 = ix2 + idot;
 742       ix4 = ix3 + idot;
 743       if (na == 0) {
 744         pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3],
 745                &wa[ix4], isign);
 746       }
 747       else {
 748         pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3],
 749                &wa[ix4], isign);
 750       }
 751       na = 1 - na;
 752     }
 753     else {
 754       if (na == 0) {
 755         pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1],
 756               &ch[1], &wa[iw], isign);
 757       }
 758       else {
 759         pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1],
 760               &c[1], &wa[iw], isign);
 761       }
 762       if (nac != 0) {
 763         na = 1 - na;
 764       }
 765     }
 766     l1 = l2;
 767     iw += (ip - 1) * idot;
 768   }
 769   if (na != 0) {
 770     n2 = *n + *n;
 771     i_1 = n2;
 772     for (i = 1; i <= i_1; ++i) {
 773       c[i] = ch[i];
 774     }
 775   }
 776   return 0;
 777 }
 778
 779
 780 /*
 781   Public Routine
 782 */
 783
 784 /*
 785  * cfft computes the forward or backward complex discrete fourier transform.
 786  *
 787  * Input Parameters:
 788  *
 789  * n      The length of the complex sequence c. The method is
 790  *        more efficient when n is the product of small primes.
 791  *
 792  * c      A complex array of length n which contains the sequence
 793  *
 794  * wsave  a real work array which must be dimensioned at least 4n+15
 795  *        in the program that calls cfft.
 796  * isign  1 for transform, -1 for inverse transform.
 797  *        A call of cfft with isign = 1 followed by a call of cfft with
 798  *        isign = -1 will multiply the sequence by n.
 799  *
 800  * Output Parameters:
 801  *
 802  * c      For j=1,...,n
 803  *
 804  *             c(j)=the sum from k=1,...,n of
 805  *
 806  *                   c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n)
 807  *
 808  *                         where i=sqrt(-1)
 809  */
 810
 811 int cfft(int n, double *c, double *wsave, int isign)
 812 {
 813   int iw1, iw2;
 814
 815   /* Parameter adjustments */
 816   --c;
 817   --wsave;
 818
 819   /* Function Body */
 820   if (n != 1) {
 821     iw1 = n + n + 1;
 822     iw2 = iw1 + n + n;
 823     cffti1_(&n, &wsave[iw1], &wsave[iw2]);
 824     cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign);
 825   }
 826   return 0;
 827 }