/* from fftpkg package in cmlib and netlib -- translated by f2c and modified */

#include "xlisp.h"
#include "xlstat.h"

/* forward declarations */
static int cfft1_ P6H(int *, double *, double *, double *, int *, int *);
static int cffti1_ P3H(int *, double *, int *);
static int pass_ P12H(int *, int *, int *, int *, int *,
		              double *, double *, double *, double *,
		              double *, double *, int *);
static int pass2_ P6H(int *, int *, double *, double *, double *, int *);
static int pass3_ P7H(int *, int *, double *, double *, double *, double *, int *);
static int pass4_ P8H(int *, int *, double *, double *, double *, double *, double *, int *);
static int pass5_ P9H(int *, int *, double *, double *,
                      double *, double *, double *, double *, int *);


/*
  Public Routine
*/

/*
 * cfft computes the forward or backward complex discrete fourier transform. 
 *
 * Input Parameters:
 *
 * n      The length of the complex sequence c. The method is
 *        more efficient when n is the product of small primes.
 *
 * c      A complex array of length n which contains the sequence
 *
 * wsave  a real work array which must be dimensioned at least 4n+15
 *        in the program that calls cfft.
 * isign  1 for transform, -1 for inverse transform.
 *        A call of cfft with isign = 1 followed by a call of cfft with 
 *        isign = -1 will multiply the sequence by n.
 *
 * Output Parameters:
 *
 * c      For j=1,...,n
 *
 *             c(j)=the sum from k=1,...,n of
 *
 *                   c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n)
 *
 *                         where i=sqrt(-1)
 */

int cfft P4C(int, n, double *, c, double *, wsave, int, isign)
{
  int iw1, iw2;
  
  /* Parameter adjustments */
  --c;
  --wsave;

  /* Function Body */
  if (n != 1) {
    iw1 = n + n + 1;
    iw2 = iw1 + n + n;
    cffti1_(&n, &wsave[iw1], (int *) &wsave[iw2]);
    cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], (int *) &wsave[iw2], &isign);
  }
  return 0;
}

/*
  Internal Routines
*/

static int cfft1_ P6C(int *, n, double *, c, double *, ch, double *, wa, int *, ifac, int *, isign)
{
  /* System generated locals */
  int i_1;

  /* Local variables */
  int idot;
  int i, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;

  /* Parameter adjustments */
  --c;
  --ch;
  --wa;
  --ifac;
  
  /* Function Body */
  nf = ifac[2];
  na = 0;
  l1 = 1;
  iw = 1;
  i_1 = nf;
  for (k1 = 1; k1 <= i_1; ++k1) {
    ip = ifac[k1 + 2];
    l2 = ip * l1;
    ido = *n / l2;
    idot = ido + ido;
    idl1 = idot * l1;
    if (ip == 4) {
      ix2 = iw + idot;
      ix3 = ix2 + idot;
      if (na == 0) {
	pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
      }
      else {
	pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
      }
      na = 1 - na;
    }
    else if (ip == 2) {
      if (na == 0) {
	pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign);
      }
      else {
	pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign);
      }
      na = 1 - na;
    }
    else if (ip == 3) {
      ix2 = iw + idot;
      if (na == 0) {
	pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign);
      }
      else {
	pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign);
      }
      na = 1 - na;
    }
    else if (ip == 5) {
      ix2 = iw + idot;
      ix3 = ix2 + idot;
      ix4 = ix3 + idot;
      if (na == 0) {
	pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], 
	       &wa[ix4], isign);
      }
      else {
	pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], 
	       &wa[ix4], isign);
      }
      na = 1 - na;
    }
    else {
      if (na == 0) {
	pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], 
	      &ch[1], &wa[iw], isign);
      }
      else {
	pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], 
	      &c[1], &wa[iw], isign);
      }
      if (nac != 0) {
	na = 1 - na;
      }
    }
    l1 = l2;
    iw += (ip - 1) * idot;
  }
  if (na != 0) {
    n2 = *n + *n;
    i_1 = n2;
    for (i = 1; i <= i_1; ++i) {
      c[i] = ch[i];
    }
  }
  return 0;
}

static int cffti1_ P3C(int *, n, double *, wa, int *, ifac)
{
  /* Initialized data */
  static int ntryh[4] = { 3,4,2,5 };

  /* System generated locals */
  int i_1, i_2, i_3;

  /* Local variables */
  double argh;
  int idot, ntry = 0, i, j;
  double argld;
  int i1, k1, l1, l2, ib;
  double fi;
  int ld, ii, nf, ip, nl, nq, nr;
  double arg;
  int ido, ipm;
  double tpi;

  /* Parameter adjustments */
  --wa;
  --ifac;

  /* Function Body */
  nl = *n;
  nf = 0;
  j = 0;

 L101:
  ++j;
  if (j - 4 <= 0) ntry = ntryh[j - 1];
  else ntry += 2;
 L104:
  nq = nl / ntry;
  nr = nl - ntry * nq;
  if (nr != 0) goto L101;
  ++nf;
  ifac[nf + 2] = ntry;
  nl = nq;
  if (ntry == 2 && nf != 1) {
    i_1 = nf;
    for (i = 2; i <= i_1; ++i) {
      ib = nf - i + 2;
      ifac[ib + 2] = ifac[ib + 1];
    }
    ifac[3] = 2;
  }
  if (nl != 1) goto L104;

  ifac[1] = *n;
  ifac[2] = nf;
  tpi = 6.28318530717959;
  argh = tpi / (double) (*n);
  i = 2;
  l1 = 1;
  i_1 = nf;
  for (k1 = 1; k1 <= i_1; ++k1) {
    ip = ifac[k1 + 2];
    ld = 0;
    l2 = l1 * ip;
    ido = *n / l2;
    idot = ido + ido + 2;
    ipm = ip - 1;
    i_2 = ipm;
    for (j = 1; j <= i_2; ++j) {
      i1 = i;
      wa[i - 1] = 1.0;
      wa[i] = 0.0;
      ld += l1;
      fi = 0.0;
      argld = (double) ld * argh;
      i_3 = idot;
      for (ii = 4; ii <= i_3; ii += 2) {
	i += 2;
	fi += 1.0;
	arg = fi * argld;
	wa[i - 1] = cos(arg);
	wa[i] = sin(arg);
      }
      if (ip > 5) {
	wa[i1 - 1] = wa[i - 1];
	wa[i1] = wa[i];
      }
    }
    l1 = l2;
  }
  return 0;
}

static int pass_ P12C(int *, nac, int *, ido, int *, ip, int *, l1, int *, idl1,
                      double *, cc, double *, c1, double *, c2, double *, ch, double *, ch2,
                      double *, wa, int *, isign)
{
  /* System generated locals */
  int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
      c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, 
      i_1, i_2, i_3;

  /* Local variables */
  int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp;
  double wai, war;
  int ipp2;

  /* Parameter adjustments */
  cc_dim1 = *ido;
  cc_dim2 = *ip;
  cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
  cc -= cc_offset;
  c1_dim1 = *ido;
  c1_dim2 = *l1;
  c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
  c1 -= c1_offset;
  c2_dim1 = *idl1;
  c2_offset = c2_dim1 + 1;
  c2 -= c2_offset;
  ch_dim1 = *ido;
  ch_dim2 = *l1;
  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
  ch -= ch_offset;
  ch2_dim1 = *idl1;
  ch2_offset = ch2_dim1 + 1;
  ch2 -= ch2_offset;
  --wa;
  
  /* Function Body */
  idot = *ido / 2;
  ipp2 = *ip + 2;
  ipph = (*ip + 1) / 2;
  idp = *ip * *ido;
  
  if (*ido >= *l1) {
    i_1 = ipph;
    for (j = 2; j <= i_1; ++j) {
      jc = ipp2 - j;
      i_2 = *l1;
      for (k = 1; k <= i_2; ++k) {
	i_3 = *ido;
	for (i = 1; i <= i_3; ++i) {
	  ch[i + (k + j * ch_dim2) * ch_dim1] =
	    cc[i + (j + k * cc_dim2) * cc_dim1] 
	      + cc[i + (jc + k * cc_dim2) * cc_dim1];
	  ch[i + (k + jc * ch_dim2) * ch_dim1] =
	    cc[i + (j + k * cc_dim2) * cc_dim1]
	      - cc[i + (jc + k * cc_dim2) * cc_dim1];
	}
      }
    }
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      i_2 = *ido;
      for (i = 1; i <= i_2; ++i) {
	ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
      }
    }
  }
  else {
    i_1 = ipph;
    for (j = 2; j <= i_1; ++j) {
      jc = ipp2 - j;
      i_2 = *ido;
      for (i = 1; i <= i_2; ++i) {
	i_3 = *l1;
	for (k = 1; k <= i_3; ++k) {
	  ch[i + (k + j * ch_dim2) * ch_dim1] = 
	    cc[i + (j + k * cc_dim2) * cc_dim1] 
	      + cc[i + (jc + k * cc_dim2) * cc_dim1];
	  ch[i + (k + jc * ch_dim2) * ch_dim1] =
	    cc[i + (j + k * cc_dim2) * cc_dim1] 
	      - cc[i + (jc + k * cc_dim2) * cc_dim1];
	}
      }
    }
    i_1 = *ido;
    for (i = 1; i <= i_1; ++i) {
      i_2 = *l1;
      for (k = 1; k <= i_2; ++k) {
	ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
      }
    }
  }
  idl = 2 - *ido;
  inc = 0;
  i_1 = ipph;
  for (l = 2; l <= i_1; ++l) {
    lc = ipp2 - l;
    idl += *ido;
    i_2 = *idl1;
    for (ik = 1; ik <= i_2; ++ik) {
      c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] 
	* ch2[ik + (ch2_dim1 << 1)];
      c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1];
    }
    idlj = idl;
    inc += *ido;
    i_2 = ipph;
    for (j = 3; j <= i_2; ++j) {
      jc = ipp2 - j;
      idlj += inc;
      if (idlj > idp) {
	idlj -= idp;
      }
      war = wa[idlj - 1];
      wai = wa[idlj];
      i_3 = *idl1;
      for (ik = 1; ik <= i_3; ++ik) {
	c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
	c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1];
      }
    }
  }
  i_1 = ipph;
  for (j = 2; j <= i_1; ++j) {
    i_2 = *idl1;
    for (ik = 1; ik <= i_2; ++ik) {
      ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
    }
  }
  i_1 = ipph;
  for (j = 2; j <= i_1; ++j) {
    jc = ipp2 - j;
    i_2 = *idl1;
    for (ik = 2; ik <= i_2; ik += 2) {
      ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
	- c2[ik + jc * c2_dim1];
      ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] 
	+ c2[ik + jc * c2_dim1];
      ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1]
	+ c2[ik - 1 + jc * c2_dim1];
      ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1]
	- c2[ik - 1 + jc * c2_dim1];
    }
  }
  *nac = 1;
  if (*ido != 2) {
    *nac = 0;
    i_1 = *idl1;
    for (ik = 1; ik <= i_1; ++ik) {
      c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
    }
    i_1 = *ip;
    for (j = 2; j <= i_1; ++j) {
      i_2 = *l1;
      for (k = 1; k <= i_2; ++k) {
	c1[(k + j * c1_dim2) * c1_dim1 + 1] = 
	  ch[(k + j * ch_dim2) * ch_dim1 + 1];
	c1[(k + j * c1_dim2) * c1_dim1 + 2] =
	  ch[(k + j * ch_dim2) * ch_dim1 + 2];
      }
    }
    if (idot <= *l1) {
      idij = 0;
      i_1 = *ip;
      for (j = 2; j <= i_1; ++j) {
	idij += 2;
	i_2 = *ido;
	for (i = 4; i <= i_2; i += 2) {
	  idij += 2;
	  i_3 = *l1;
	  for (k = 1; k <= i_3; ++k) {
	    c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] 
	      * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij]
		* ch[i + (k + j * ch_dim2) * ch_dim1];
	    c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] 
	      * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] 
		* ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
	  }
	}
      }
    } 
    else {
      idj = 2 - *ido;
      i_1 = *ip;
      for (j = 2; j <= i_1; ++j) {
	idj += *ido;
	i_2 = *l1;
	for (k = 1; k <= i_2; ++k) {
	  idij = idj;
	  i_3 = *ido;
	  for (i = 4; i <= i_3; i += 2) {
	    idij += 2;
	    c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] 
	      * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] 
		+ *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1];
	    c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] 
	      * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] 
		* ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
	  }
	}
      }
    }
  }
  return 0;
}

static int pass2_ P6C(int *, ido, int *, l1, double *, cc, double *, ch, double *, wa1, int *, isign)
{
  /* System generated locals */
  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
  
  /* Local variables */
  int i, k;
  double ti2, tr2;
  
  /* Parameter adjustments */
  cc_dim1 = *ido;
  cc_offset = cc_dim1 * 3 + 1;
  cc -= cc_offset;
  ch_dim1 = *ido;
  ch_dim2 = *l1;
  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
  ch -= ch_offset;
  --wa1;
  
  /* Function Body */
  if (*ido <= 2) {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
	+ cc[((k << 1) + 2) * cc_dim1 + 1];
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
	- cc[((k << 1) + 2) * cc_dim1 + 1];
      ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
	+ cc[((k << 1) + 2) * cc_dim1 + 2];
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
	- cc[((k << 1) + 2) * cc_dim1 + 2];
    }
  }
  else {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      i_2 = *ido;
      for (i = 2; i <= i_2; i += 2) {
	ch[i - 1 + (k + ch_dim2) * ch_dim1] 
	  = cc[i - 1 + ((k << 1) + 1) * cc_dim1] 
	    + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
	tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] 
	  - cc[i - 1 + ((k << 1) + 2) * cc_dim1];
	ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1] 
	  + cc[i + ((k << 1) + 2) * cc_dim1];
	ti2 = cc[i + ((k << 1) + 1) * cc_dim1] 
	  - cc[i + ((k << 1) + 2) * cc_dim1];
	ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 
	  - *isign * wa1[i] * tr2;
	ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 
	  + *isign * wa1[i] * ti2;
      }
    }
  }
  return 0;
}

static int pass3_ P7C(int *, ido, int *, l1,
                      double *, cc, double *, ch, double *, wa1, double *, wa2,
                      int *, isign)
{
  /* System generated locals */
  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;

  /* Local variables */
  double taui, taur;
  int i, k;
  double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;

  /* Parameter adjustments */
  cc_dim1 = *ido;
  cc_offset = (cc_dim1 << 2) + 1;
  cc -= cc_offset;
  ch_dim1 = *ido;
  ch_dim2 = *l1;
  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
  ch -= ch_offset;
  --wa1;
  --wa2;

  /* Function Body */
  taur = -.5;
  taui = -(*isign) * .866025403784439;
  if (*ido == 2) {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
      cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
      ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
      ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
      ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
      ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
      cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] 
		    - cc[(k * 3 + 3) * cc_dim1 + 1]);
      ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] 
		    - cc[(k * 3 + 3) * cc_dim1 + 2]);
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
    }
  }
  else {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      i_2 = *ido;
      for (i = 2; i <= i_2; i += 2) {
	tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] 
	  + cc[i - 1 + (k * 3 + 3) * cc_dim1];
	cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
	ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1)
						 * cc_dim1] + tr2;
	ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1];
	ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
	ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2;
	cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] 
		      - cc[i - 1 + (k * 3 + 3) * cc_dim1]);
	ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] 
		      - cc[i + (k * 3 + 3) * cc_dim1]);
	dr2 = cr2 - ci3;
	dr3 = cr2 + ci3;
	di2 = ci2 + cr3;
	di3 = ci2 - cr3;
	ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 
	  - *isign * wa1[i] * dr2;
	ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 
	  + *isign * wa1[i] * di2;
	ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 
	  - *isign * wa2[i] * dr3;
	ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 
	  + *isign * wa2[i] * di3;
      }
    }
  }
  return 0;
}

static int pass4_ P8C(int *, ido, int *, l1,
                      double *, cc, double *, ch, double *, wa1, double *, wa2, double *, wa3,
                      int *, isign)
{
  /* System generated locals */
  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;

  /* Local variables */
  int i, k;
  double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;

  /* Parameter adjustments */
  cc_dim1 = *ido;
  cc_offset = cc_dim1 * 5 + 1;
  cc -= cc_offset;
  ch_dim1 = *ido;
  ch_dim2 = *l1;
  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
  ch -= ch_offset;
  --wa1;
  --wa2;
  --wa3;
  
  /* Function Body */
  if (*ido == 2) {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] 
	- cc[((k << 2) + 3) * cc_dim1 + 2];
      ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] 
	+ cc[((k << 2) + 3) * cc_dim1 + 2];
      tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2] 
		      - cc[((k << 2) + 4) * cc_dim1 + 2]);
      ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] 
	+ cc[((k << 2) + 4) * cc_dim1 + 2];
      tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] 
	- cc[((k << 2) + 3) * cc_dim1 + 1];
      tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] 
	+ cc[((k << 2) + 3) * cc_dim1 + 1];
      ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1] 
		      - cc[((k << 2) + 2) * cc_dim1 + 1]);
      tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] 
	+ cc[((k << 2) + 4) * cc_dim1 + 1];
      ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
      ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
      ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
      ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
    }
  }
  else {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      i_2 = *ido;
      for (i = 2; i <= i_2; i += 2) {
	ti1 = cc[i + ((k << 2) + 1) * cc_dim1] 
	  - cc[i + ((k << 2) + 3) * cc_dim1];
	ti2 = cc[i + ((k << 2) + 1) * cc_dim1] 
	  + cc[i + ((k << 2) + 3) * cc_dim1];
	ti3 = cc[i + ((k << 2) + 2) * cc_dim1] 
	  + cc[i + ((k << 2) + 4) * cc_dim1];
	tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1] 
			- cc[i + ((k << 2) + 4) * cc_dim1]);
	tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] 
	  - cc[i - 1 + ((k << 2) + 3) * cc_dim1];
	tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] 
	  + cc[i - 1 + ((k << 2) + 3) * cc_dim1];
	ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1]
			- cc[i - 1 + ((k << 2) + 2) * cc_dim1]);
	tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] 
	  + cc[i - 1 + ((k << 2) + 4) * cc_dim1];
	ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
	cr3 = tr2 - tr3;
	ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
	ci3 = ti2 - ti3;
	cr2 = tr1 + tr4;
	cr4 = tr1 - tr4;
	ci2 = ti1 + ti4;
	ci4 = ti1 - ti4;
	ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 
	  + *isign * wa1[i] * ci2;
	ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 
	  - *isign * wa1[i] * cr2;
	ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3
	  + *isign * wa2[i] * ci3;
	ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 
	  - *isign * wa2[i] * cr3;
	ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 
	  + *isign * wa3[i] * ci4;
	ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 
	  - *isign * wa3[i] * cr4;
      }
    }
  }
  return 0;
}

static int pass5_ P9C(int *, ido, int *, l1,
                      double *, cc, double *, ch,
                      double *, wa1, double *, wa2, double *, wa3, double *, wa4,
                      int *, isign)
{
  /* System generated locals */
  int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;

  /* Local variables */
  int i, k;
  double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
         ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11,
         tr12;

  /* Parameter adjustments */
  cc_dim1 = *ido;
  cc_offset = cc_dim1 * 6 + 1;
  cc -= cc_offset;
  ch_dim1 = *ido;
  ch_dim2 = *l1;
  ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
  ch -= ch_offset;
  --wa1;
  --wa2;
  --wa3;
  --wa4;
  
  /* Function Body */
  tr11 = .309016994374947;
  ti11 = -(*isign) * .951056516295154;
  tr12 = -.809016994374947;
  ti12 = -(*isign) * .587785252292473;
  if (*ido == 2) {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
      ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
      ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
      ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
      tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
      tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
      tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
      tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
      ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] 
	+ tr2 + tr3;
      ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] 
	+ ti2 + ti3;
      cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
      ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
      cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
      ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
      cr5 = ti11 * tr5 + ti12 * tr4;
      ci5 = ti11 * ti5 + ti12 * ti4;
      cr4 = ti12 * tr5 - ti11 * tr4;
      ci4 = ti12 * ti5 - ti11 * ti4;
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
      ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
      ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
      ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
      ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
      ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
      ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
    }
  }
  else {
    i_1 = *l1;
    for (k = 1; k <= i_1; ++k) {
      i_2 = *ido;
      for (i = 2; i <= i_2; i += 2) {
	ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1];
	ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1];
	ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1];
	ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1];
	tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] 
	  - cc[i - 1 + (k * 5 + 5) * cc_dim1];
	tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] 
	  + cc[i - 1 + (k * 5 + 5) * cc_dim1];
	tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] 
	  - cc[i - 1 + (k * 5 + 4) * cc_dim1];
	tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] 
	  + cc[i - 1 + (k * 5 + 4) * cc_dim1];
	ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) 
						 * cc_dim1] + tr2 + tr3;
	ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) 
					     * cc_dim1] + ti2 + ti3;
	cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
	ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;

	cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
	ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;

	cr5 = ti11 * tr5 + ti12 * tr4;
	ci5 = ti11 * ti5 + ti12 * ti4;
	cr4 = ti12 * tr5 - ti11 * tr4;
	ci4 = ti12 * ti5 - ti11 * ti4;
	dr3 = cr3 - ci4;
	dr4 = cr3 + ci4;
	di3 = ci3 + cr4;
	di4 = ci3 - cr4;
	dr5 = cr2 + ci5;
	dr2 = cr2 - ci5;
	di5 = ci2 - cr5;
	di2 = ci2 + cr5;
	ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 
	  + *isign * wa1[i] * di2;
	ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
	  - *isign * wa1[i] * dr2;
	ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 
	  + *isign * wa2[i] * di3;
	ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 
	  - *isign * wa2[i] * dr3;
	ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 
	  + *isign * wa3[i] * di4;
	ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4
	  - *isign * wa3[i] * dr4;
	ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 
	  + *isign * wa4[i] * di5;
	ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 
	  - *isign * wa4[i] * dr5;
      }
    }
  }
  return 0;
}