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