Logo Search packages:      
Sourcecode: tela version File versions  Download package

vecfftpack.c

/*
 * This file is part of tela the Tensor Language.
 * Copyright (c) 1994-1996 Pekka Janhunen
 */

/*
  vecfftpack.c : A set of vectorized FFT routines (many transforms with one call).
  Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985).
  As compared to FFTPACK routines, these functions take three extra integer parameters:
  lot, step and jump.
  Lot is the number of transforms to perform with the call.
  Step is the address difference between consecutive components in each data vector.
  Jump is the address difference between corresponding elements of consecutive
  vectors.
  These parameters are arbitrary and thus it is possible to take the transforms along
  rows or columns of a matrix.

  Assumption throughout is that sizeof(Tint) <= sizeof(Treal).
  E.g. on modern 64-bit RISC machines you cannot define Treal=float, Tint=long.
  This is checked for in the code.

  Compilation outside Tela: cc -DSTANDALONE vecfftpack.c.
  Vectorized compilation on Cray: cc -DUNICOS -DSTANDALONE vecfftpack.c.
  If you don't define STANDALONE it tries to include "deftyp.H".
  
  This file is also C++ compatible.
  
  Pekka Janhunen 23.2.1995
*/


/* isign is +1 for backward and -1 for forward transforms */

#include <math.h>
#include <stdio.h>

#ifdef STANDALONE
#  define Treal double
#  define Tint int
#else
#  include "deftyp.h"
#endif

#define ref(u,a) u[p*jump+(a)*step]
#define ref0(u,a) u[2*p*jump+2*((a)/2)*step+(a)%2]
#define ref1(u,a) u[2*p*jump+2*(((a)+1)/2)*step+((a)+1)%2]

#define MAXFAC 13

#define VECARGS_FORMAL , int lot, int step, int jump
#define VECARGS_ACTUAL ,lot,step,jump
#define VECVARS Tint p
#ifdef UNICOS
#  define VECLOOP _Pragma("_CRI ivdep"); for (p=0; p<lot; p++)
#else
#  define VECLOOP for (p=0; p<lot; p++)
#endif

static const Treal pi       = 3.14159265358979;
static const Treal twopi    = 6.28318530717959;
static const Treal taur     = -0.5;
static const Treal taui     = 0.866025403784439;            /* sqrt(3)/2 */
static const Treal tr11     = 0.309016994374947;            /* sin(pi/10) */
static const Treal ti11     = 0.951056516295154;            /* cos(pi/10) */
static const Treal tr12     = -0.809016994374947;           /* -cos(pi/5) */
static const Treal ti12     = 0.587785252292473;            /* sin(pi/5) */
static const Treal sqrt2    = 1.4142135623731;
static const Treal tsqrt2   = 2.82842712474619;             /* 2*sqrt(2) */
static const Treal invsqrt2 = 0.7071067811865475;           /* 1/sqrt(2) */
static const Treal sqrt3    = 1.73205080756888;

#ifdef __cplusplus
extern "C" {
#endif

#ifdef UNICOS
#pragma _CRI noinline passf2,passf3,passf4,passf5,passf
#pragma _CRI noinline radf2,radf3,radf4,radf5,radfg
#pragma _CRI noinline radb2,radb3,radb4,radb5,radbg
#pragma _CRI noinline cfftf1,rfftf1
#pragma _CRI inline cosqf1,cosqb1
#pragma _CRI inline factorize,rfftb1,rffti1,sint1
#endif

/* ----------------------------------------------------------------------
   passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
   ---------------------------------------------------------------------- */

static void passf2(Tint ido, Tint l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign VECARGS_FORMAL)
/* isign==+1 for backward transform */
{
      VECVARS;
    Tint i, k, ah, ac;
    Treal ti2, tr2;
    if (ido <= 2) {
            for (k=0; k<l1; k++) {
                  ah = k*ido;
                  ac = 2*k*ido;
                  VECLOOP {
                  ref0(ch,ah)           = ref0(cc,ac) + ref0(cc,ac + ido);
                  ref0(ch,ah + ido*l1)  = ref0(cc,ac) - ref0(cc,ac + ido);
                  ref1(ch,ah)           = ref1(cc,ac) + ref1(cc,ac + ido);
                  ref1(ch,ah + ido*l1)  = ref1(cc,ac) - ref1(cc,ac + ido);
                  }
            }
      } else {
            for (k=0; k<l1; k++) {
                  for (i=0; i<ido-1; i+=2) {
                        ah = i + k*ido;
                        ac = i + 2*k*ido;
                        VECLOOP {
                        ref0(ch,ah)         = ref0(cc,ac) + ref0(cc,ac + ido);
                        tr2                 = ref0(cc,ac) - ref0(cc,ac + ido);
                        ref1(ch,ah)         = ref1(cc,ac) + ref1(cc,ac + ido);
                        ti2                 = ref1(cc,ac) - ref1(cc,ac + ido);
                        ref1(ch,ah+l1*ido)  = wa1[i]*ti2 + isign*wa1[i+1]*tr2;
                        ref0(ch,ah+l1*ido)  = wa1[i]*tr2 - isign*wa1[i+1]*ti2;
                        }
                  }
            }
      }
} /* passf2 */


static void passf3(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                           const Treal wa1[], const Treal wa2[], int isign VECARGS_FORMAL)
/* isign==+1 for backward transform */
{
      VECVARS;
    Tint i, k, ac, ah;
    Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
    if (ido == 2) {
            for (k=1; k<=l1; k++) {
                  VECLOOP {
                  ac = (3*k - 2)*ido;
                  tr2 = ref0(cc,ac) + ref0(cc,ac + ido);
                  cr2 = ref0(cc,ac - ido) + taur*tr2;
                  ah = (k - 1)*ido;
                  ref0(ch, ah) = ref0(cc,ac - ido) + tr2;

                  ti2 = ref1(cc,ac) + ref1(cc,ac + ido);
                  ci2 = ref1(cc,ac - ido) + taur*ti2;
                  ref1(ch, ah) = ref1(cc,ac - ido) + ti2;

                  cr3 = isign*taui*(ref0(cc,ac) - ref0(cc,ac + ido));
                  ci3 = isign*taui*(ref1(cc,ac) - ref1(cc,ac + ido));
                  ref0(ch, ah + l1*ido) = cr2 - ci3;
                  ref0(ch, ah + 2*l1*ido) = cr2 + ci3;
                  ref1(ch, ah + l1*ido) = ci2 + cr3;
                  ref1(ch, ah + 2*l1*ido) = ci2 - cr3;
                  }
            }
      } else {
            for (k=1; k<=l1; k++) {
                  for (i=0; i<ido-1; i+=2) {
                        VECLOOP {
                        ac = i + (3*k - 2)*ido;
                        tr2 = ref0(cc,ac) + ref0(cc,ac + ido);
                        cr2 = ref0(cc,ac - ido) + taur*tr2;
                        ah = i + (k-1)*ido;
                        ref0(ch, ah) = ref0(cc,ac - ido) + tr2;
                        ti2 = ref1(cc,ac) + ref1(cc,ac + ido);
                        ci2 = ref1(cc,ac - ido) + taur*ti2;
                        ref1(ch, ah) = ref1(cc,ac - ido) + ti2;
                        cr3 = isign*taui*(ref0(cc,ac) - ref0(cc,ac + ido));
                        ci3 = isign*taui*(ref1(cc,ac) - ref1(cc,ac + ido));
                        dr2 = cr2 - ci3;
                        dr3 = cr2 + ci3;
                        di2 = ci2 + cr3;
                        di3 = ci2 - cr3;
                        ref1(ch, ah + l1*ido) = wa1[i]*di2 + isign*wa1[i+1]*dr2;
                        ref0(ch, ah + l1*ido) = wa1[i]*dr2 - isign*wa1[i+1]*di2;
                        ref1(ch, ah + 2*l1*ido) = wa2[i]*di3 + isign*wa2[i+1]*dr3;
                        ref0(ch, ah + 2*l1*ido) = wa2[i]*dr3 - isign*wa2[i+1]*di3;
                        }
                  }
            }
      }
} /* passf3 */


static void passf4(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                           const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign VECARGS_FORMAL)
/* isign == -1 for forward transform and +1 for backward transform */
{
      VECVARS;
    Tint i, k, ac, ah;
    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
    if (ido == 2) {
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ac = 4*k*ido;
                  ti1 = ref1(cc,ac) - ref1(cc,ac + 2*ido);
                  ti2 = ref1(cc,ac) + ref1(cc,ac + 2*ido);
                  tr4 = ref1(cc,ac + 3*ido) - ref1(cc,ac + ido);
                  ti3 = ref1(cc,ac + ido) + ref1(cc,ac + 3*ido);
                  tr1 = ref0(cc,ac) - ref0(cc,ac + 2*ido);
                  tr2 = ref0(cc,ac) + ref0(cc,ac + 2*ido);
                  ti4 = ref0(cc,ac + ido) - ref0(cc,ac + 3*ido);
                  tr3 = ref0(cc,ac + ido) + ref0(cc,ac + 3*ido);
                  ah = k*ido;
                  ref0(ch, ah) = tr2 + tr3;
                  ref0(ch, ah + 2*l1*ido) = tr2 - tr3;
                  ref1(ch, ah) = ti2 + ti3;
                  ref1(ch, ah + 2*l1*ido) = ti2 - ti3;
                  ref0(ch, ah + l1*ido) = tr1 + isign*tr4;
                  ref0(ch, ah + 3*l1*ido) = tr1 - isign*tr4;
                  ref1(ch, ah + l1*ido) = ti1 + isign*ti4;
                  ref1(ch, ah + 3*l1*ido) = ti1 - isign*ti4;
                  }
            }
      } else {
            for (k=0; k<l1; k++) {
                  for (i=0; i<ido-1; i+=2) {
                        VECLOOP {
                        ac = i + 1 + 4*k*ido;
                        ti1 = ref0(cc,ac) - ref0(cc,ac + 2*ido);
                        ti2 = ref0(cc,ac) + ref0(cc,ac + 2*ido);
                        ti3 = ref0(cc,ac + ido) + ref0(cc,ac + 3*ido);
                        tr4 = ref0(cc,ac + 3*ido) - ref0(cc,ac + ido);
                        tr1 = ref0(cc,ac - 1) - ref0(cc,ac + 2*ido - 1);
                        tr2 = ref0(cc,ac - 1) + ref0(cc,ac + 2*ido - 1);
                        ti4 = ref0(cc,ac + ido - 1) - ref0(cc,ac + 3*ido - 1);
                        tr3 = ref0(cc,ac + ido - 1) + ref0(cc,ac + 3*ido - 1);
                        ah = i + k*ido;
                        ref0(ch, ah) = tr2 + tr3;
                        cr3 = tr2 - tr3;
                        ref0(ch, ah + 1) = ti2 + ti3;
                        ci3 = ti2 - ti3;
                        cr2 = tr1 + isign*tr4;
                        cr4 = tr1 - isign*tr4;
                        ci2 = ti1 + isign*ti4;
                        ci4 = ti1 - isign*ti4;
                        ref0(ch, ah + l1*ido) = wa1[i]*cr2 - isign*wa1[i + 1]*ci2;
                        ref0(ch, ah + l1*ido + 1) = wa1[i]*ci2 + isign*wa1[i + 1]*cr2;
                        ref0(ch, ah + 2*l1*ido) = wa2[i]*cr3 - isign*wa2[i + 1]*ci3;
                        ref0(ch, ah + 2*l1*ido + 1) = wa2[i]*ci3 + isign*wa2[i + 1]*cr3;
                        ref0(ch, ah + 3*l1*ido) = wa3[i]*cr4 -isign*wa3[i + 1]*ci4;
                        ref0(ch, ah + 3*l1*ido + 1) = wa3[i]*ci4 + isign*wa3[i + 1]*cr4;
                        }
                  }
            }
      }
} /* passf4 */


static void passf5(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                           const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[],
                           int isign VECARGS_FORMAL)
/* isign == -1 for forward transform and +1 for backward transform */
{
      VECVARS;
    Tint i, k, ac, ah;
    Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
           ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
    if (ido == 2) {
            for (k = 1; k <= l1; ++k) {
                  VECLOOP {
                  ac = (5*k - 4)*ido + 1;
                  ti5 = ref0(cc,ac) - ref0(cc,ac + 3*ido);
                  ti2 = ref0(cc,ac) + ref0(cc,ac + 3*ido);
                  ti4 = ref0(cc,ac + ido) - ref0(cc,ac + 2*ido);
                  ti3 = ref0(cc,ac + ido) + ref0(cc,ac + 2*ido);
                  tr5 = ref0(cc,ac - 1) - ref0(cc,ac + 3*ido - 1);
                  tr2 = ref0(cc,ac - 1) + ref0(cc,ac + 3*ido - 1);
                  tr4 = ref0(cc,ac + ido - 1) - ref0(cc,ac + 2*ido - 1);
                  tr3 = ref0(cc,ac + ido - 1) + ref0(cc,ac + 2*ido - 1);
                  ah = (k - 1)*ido;
                  ref0(ch, ah) = ref0(cc,ac - ido - 1) + tr2 + tr3;
                  ref0(ch, ah + 1) = ref0(cc,ac - ido) + ti2 + ti3;
                  cr2 = ref0(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;
                  ci2 = ref0(cc,ac - ido) + tr11*ti2 + tr12*ti3;
                  cr3 = ref0(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;
                  ci3 = ref0(cc,ac - ido) + tr12*ti2 + tr11*ti3;
                  cr5 = isign*(ti11*tr5 + ti12*tr4);
                  ci5 = isign*(ti11*ti5 + ti12*ti4);
                  cr4 = isign*(ti12*tr5 - ti11*tr4);
                  ci4 = isign*(ti12*ti5 - ti11*ti4);
                  ref0(ch, ah + l1*ido) = cr2 - ci5;
                  ref0(ch, ah + 4*l1*ido) = cr2 + ci5;
                  ref0(ch, ah + l1*ido + 1) = ci2 + cr5;
                  ref0(ch, ah + 2*l1*ido + 1) = ci3 + cr4;
                  ref0(ch, ah + 2*l1*ido) = cr3 - ci4;
                  ref0(ch, ah + 3*l1*ido) = cr3 + ci4;
                  ref0(ch, ah + 3*l1*ido + 1) = ci3 - cr4;
                  ref0(ch, ah + 4*l1*ido + 1) = ci2 - cr5;
                  }
            }
      } else {
            for (k=1; k<=l1; k++) {
                  for (i=0; i<ido-1; i+=2) {
                        VECLOOP {
                        ac = i + 1 + (k*5 - 4)*ido;
                        ti5 = ref0(cc,ac) - ref0(cc,ac + 3*ido);
                        ti2 = ref0(cc,ac) + ref0(cc,ac + 3*ido);
                        ti4 = ref0(cc,ac + ido) - ref0(cc,ac + 2*ido);
                        ti3 = ref0(cc,ac + ido) + ref0(cc,ac + 2*ido);
                        tr5 = ref0(cc,ac - 1) - ref0(cc,ac + 3*ido - 1);
                        tr2 = ref0(cc,ac - 1) + ref0(cc,ac + 3*ido - 1);
                        tr4 = ref0(cc,ac + ido - 1) - ref0(cc,ac + 2*ido - 1);
                        tr3 = ref0(cc,ac + ido - 1) + ref0(cc,ac + 2*ido - 1);
                        ah = i + (k - 1)*ido;
                        ref0(ch, ah) = ref0(cc,ac - ido - 1) + tr2 + tr3;
                        ref0(ch, ah + 1) = ref0(cc,ac - ido) + ti2 + ti3;
                        cr2 = ref0(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3;

                        ci2 = ref0(cc,ac - ido) + tr11*ti2 + tr12*ti3;
                        cr3 = ref0(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3;

                        ci3 = ref0(cc,ac - ido) + tr12*ti2 + tr11*ti3;
                        cr5 = isign*(ti11*tr5 + ti12*tr4);
                        ci5 = isign*(ti11*ti5 + ti12*ti4);
                        cr4 = isign*(ti12*tr5 - ti11*tr4);
                        ci4 = isign*(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;
                        ref0(ch, ah + l1*ido) = wa1[i]*dr2 - isign*wa1[i+1]*di2;
                        ref0(ch, ah + l1*ido + 1) = wa1[i]*di2 + isign*wa1[i+1]*dr2;
                        ref0(ch, ah + 2*l1*ido) = wa2[i]*dr3 - isign*wa2[i+1]*di3;
                        ref0(ch, ah + 2*l1*ido + 1) = wa2[i]*di3 + isign*wa2[i+1]*dr3;
                        ref0(ch, ah + 3*l1*ido) = wa3[i]*dr4 - isign*wa3[i+1]*di4;
                        ref0(ch, ah + 3*l1*ido + 1) = wa3[i]*di4 + isign*wa3[i+1]*dr4;
                        ref0(ch, ah + 4*l1*ido) = wa4[i]*dr5 - isign*wa4[i+1]*di5;
                        ref0(ch, ah + 4*l1*ido + 1) = wa4[i]*di5 + isign*wa4[i+1]*dr5;
                        }
                  }
            }
      }
} /* passf5 */


static void passf(Tint *nac, Tint ido, Tint ip, Tint l1, Tint idl1,
                          Treal cc[], Treal ch[],
                          const Treal wa[], int isign VECARGS_FORMAL)
/* isign is -1 for forward transform and +1 for backward transform */
{
      VECVARS;
    Tint idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
    Treal wai, war;

    idot = ido / 2;
    nt = ip*idl1;
    ipph = (ip + 1) / 2;
    idp = ip*ido;
    if (ido >= l1) {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  for (k=0; k<l1; k++) {
                        for (i=0; i<ido; i++) {
                              VECLOOP {
                              ref0(ch, i + (k + j*l1)*ido) =
                                    ref0(cc,i + (j + k*ip)*ido) + ref0(cc,i + (jc + k*ip)*ido);
                              ref0(ch, i + (k + jc*l1)*ido) =
                                    ref0(cc,i + (j + k*ip)*ido) - ref0(cc,i + (jc + k*ip)*ido);
                              }
                        }
                  }
            }
            for (k=0; k<l1; k++)
                  for (i=0; i<ido; i++) {
                        VECLOOP
                              ref0(ch, i + k*ido) = ref0(cc,i + k*ip*ido);
                  }
      } else {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  for (i=0; i<ido; i++) {
                        for (k=0; k<l1; k++) {
                              VECLOOP {
                              ref0(ch, i + (k + j*l1)*ido) = ref0(cc,i + (j + k*ip)*ido) + ref0(cc,i + (jc + k*ip)*ido);
                              ref0(ch, i + (k + jc*l1)*ido) = ref0(cc,i + (j + k*ip)*ido) - ref0(cc,i + (jc + k*ip)*ido);
                              }
                        }
                  }
            }
            for (i=0; i<ido; i++)
                  for (k=0; k<l1; k++) {
                        VECLOOP
                              ref0(ch, i + k*ido) = ref0(cc,i + k*ip*ido);
                  }
      }
      
      idl = 2 - ido;
    inc = 0;
    for (l=1; l<ipph; l++) {
            lc = ip - l;
            idl += ido;
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP {
                  ref0(cc, ik + l*idl1) = ref0(ch, ik) + wa[idl - 2]*ref0(ch, ik + idl1);
                  ref0(cc, ik + lc*idl1) = isign*wa[idl-1]*ref0(ch, ik + (ip-1)*idl1);
                  }
            }
            idlj = idl;
            inc += ido;
            for (j=2; j<ipph; j++) {
                  jc = ip - j;
                  idlj += inc;
                  if (idlj > idp) idlj -= idp;
                  war = wa[idlj - 2];
                  wai = wa[idlj-1];
                  for (ik=0; ik<idl1; ik++) {
                  VECLOOP {
                        ref0(cc, ik + l*idl1) += war*ref0(ch, ik + j*idl1);
                        ref0(cc, ik + lc*idl1) += isign*wai*ref0(ch, ik + jc*idl1);
                  }
                  }
            }
    }
    for (j=1; j<ipph; j++)
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP
                        ref0(ch, ik) += ref0(ch, ik + j*idl1);
            }
    for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (ik=1; ik<idl1; ik+=2) {
                  VECLOOP {
                  ref0(ch, ik - 1 + j*idl1) = ref0(cc, ik - 1 + j*idl1) - ref0(cc, ik + jc*idl1);
                  ref0(ch, ik - 1 + jc*idl1) = ref0(cc, ik - 1 + j*idl1) + ref0(cc, ik + jc*idl1);
                  ref0(ch, ik + j*idl1) = ref0(cc, ik + j*idl1) + ref0(cc, ik - 1 + jc*idl1);
                  ref0(ch, ik + jc*idl1) = ref0(cc, ik + j*idl1) - ref0(cc, ik - 1 + jc*idl1);
                  }
            }
    }
    if (ido == 2) {*nac = 1; return;}
    *nac = 0;
    for (ik=0; ik<idl1; ik++) {
            VECLOOP
                  ref0(cc, ik) = ref0(ch, ik);
      }
    for (j=1; j<ip; j++) {
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ref0(cc, (k + j*l1)*ido + 0) = ref0(ch, (k + j*l1)*ido + 0);
                  ref0(cc, (k + j*l1)*ido + 1) = ref0(ch, (k + j*l1)*ido + 1);
                  }
            }
    }
    if (idot <= l1) {
            idij = 0;
            for (j=1; j<ip; j++) {
                  idij += 2;
                  for (i=3; i<ido; i+=2) {
                        idij += 2;
                        for (k=0; k<l1; k++) {
                              VECLOOP {
                              ref0(cc, i - 1 + (k + j*l1)*ido) =
                                    wa[idij - 2]*ref0(ch, i - 1 + (k + j*l1)*ido) -
                                    isign*wa[idij-1]*ref0(ch, i + (k + j*l1)*ido);
                              ref0(cc, i + (k + j*l1)*ido) =
                                    wa[idij - 2]*ref0(ch, i + (k + j*l1)*ido) +
                                    isign*wa[idij-1]*ref0(ch, i - 1 + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      } else {
            idj = 2 - ido;
            for (j=1; j<ip; j++) {
                  idj += ido;
                  for (k = 0; k < l1; k++) {
                        idij = idj;
                        for (i=3; i<ido; i+=2) {
                              idij += 2;
                              VECLOOP {
                              ref0(cc, i - 1 + (k + j*l1)*ido) =
                                    wa[idij - 2]*ref0(ch, i - 1 + (k + j*l1)*ido) -
                                    isign*wa[idij-1]*ref0(ch, i + (k + j*l1)*ido);
                              ref0(cc, i + (k + j*l1)*ido) =
                                    wa[idij - 2]*ref0(ch, i + (k + j*l1)*ido) +
                                    isign*wa[idij-1]*ref0(ch, i - 1 + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      }
} /* passf */


/* ----------------------------------------------------------------------
   radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
   Real FFT passes fwd and bwd.
   ---------------------------------------------------------------------- */

static void radf2(Tint ido, Tint l1, const Treal cc[], Treal ch[], const Treal wa1[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ti2, tr2;
    for (k=0; k<l1; k++) {
            VECLOOP {
            ref(ch, 2*k*ido) =
                  ref(cc,k*ido) + ref(cc,(k + l1)*ido);
            ref(ch, (2*k+1)*ido + ido-1) =
                  ref(cc,k*ido) - ref(cc,(k + l1)*ido);
            }
    }
    if (ido < 2) return;
      if (ido != 2) {
            for (k=0; k<l1; k++) {
                  for (i=2; i<ido; i+=2) {
                        ic = ido - i;
                        VECLOOP {
                        tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido);
                        ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido);
                        ref(ch, i + 2*k*ido) = ref(cc,i + k*ido) + ti2;
                        ref(ch, ic + (2*k+1)*ido) = ti2 - ref(cc,i + k*ido);
                        ref(ch, i - 1 + 2*k*ido) = ref(cc,i - 1 + k*ido) + tr2;
                        ref(ch, ic - 1 + (2*k+1)*ido) = ref(cc,i - 1 + k*ido) - tr2;
                        }
                  }
            }
            if (ido % 2 == 1) return;
      }
      for (k=0; k<l1; k++) {
            VECLOOP {
            ref(ch, (2*k+1)*ido) = -ref(cc,ido-1 + (k + l1)*ido);
            ref(ch, ido-1 + 2*k*ido) = ref(cc,ido-1 + k*ido);
            }
      }
} /* radf2 */


static void radb2(Tint ido, Tint l1, const Treal cc[], Treal ch[], const Treal wa1[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ti2, tr2;
    for (k=0; k<l1; k++) {
            VECLOOP {
            ref(ch, k*ido) =
                  ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido);
            ref(ch, (k + l1)*ido) =
                  ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido);
            }
    }
    if (ido < 2) return;
    if (ido != 2) {
            for (k = 0; k < l1; ++k) {
                  for (i = 2; i < ido; i += 2) {
                        ic = ido - i;
                        VECLOOP {
                        ref(ch, i-1 + k*ido) =
                              ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido);
                        tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido);
                        ref(ch, i + k*ido) =
                              ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido);
                        ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido);
                        ref(ch, i-1 + (k + l1)*ido) =
                              wa1[i - 2]*tr2 - wa1[i - 1]*ti2;
                        ref(ch, i + (k + l1)*ido) =
                              wa1[i - 2]*ti2 + wa1[i - 1]*tr2;
                        }
                  }
            }
            if (ido % 2 == 1) return;
      }
    for (k = 0; k < l1; k++) {
            VECLOOP {
            ref(ch, ido-1 + k*ido) = 2*ref(cc,ido-1 + 2*k*ido);
            ref(ch, ido-1 + (k + l1)*ido) = -2*ref(cc,(2*k+1)*ido);
            }
    }
} /* radb2 */


static void radf3(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
    for (k=0; k<l1; k++) {
            VECLOOP {
            cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido);
            ref(ch, 3*k*ido) = ref(cc,k*ido) + cr2;
            ref(ch, (3*k+2)*ido) = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido));
            ref(ch, ido-1 + (3*k + 1)*ido) = ref(cc,k*ido) + taur*cr2;
            }
    }
    if (ido == 1) return;
    for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                  ic = ido - i;
                  VECLOOP {
                  dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + 
                        wa1[i - 1]*ref(cc,i + (k + l1)*ido);
                  di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
                  dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido);
                  di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido);
                  cr2 = dr2 + dr3;
                  ci2 = di2 + di3;
                  ref(ch, i - 1 + 3*k*ido) = ref(cc,i - 1 + k*ido) + cr2;
                  ref(ch, i + 3*k*ido) = ref(cc,i + k*ido) + ci2;
                  tr2 = ref(cc,i - 1 + k*ido) + taur*cr2;
                  ti2 = ref(cc,i + k*ido) + taur*ci2;
                  tr3 = taui*(di2 - di3);
                  ti3 = taui*(dr3 - dr2);
                  ref(ch, i - 1 + (3*k + 2)*ido) = tr2 + tr3;
                  ref(ch, ic - 1 + (3*k + 1)*ido) = tr2 - tr3;
                  ref(ch, i + (3*k + 2)*ido) = ti2 + ti3;
                  ref(ch, ic + (3*k + 1)*ido) = ti3 - ti2;
                  }
            }
    }
} /* radf3 */


static void radb3(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
    for (k=0; k<l1; k++) {
            VECLOOP {
            tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido);
            cr2 = ref(cc,3*k*ido) + taur*tr2;
            ref(ch, k*ido) = ref(cc,3*k*ido) + tr2;
            ci3 = 2*taui*ref(cc,(3*k + 2)*ido);
            ref(ch, (k + l1)*ido) = cr2 - ci3;
            ref(ch, (k + 2*l1)*ido) = cr2 + ci3;
            }
    }
    if (ido == 1) return;
    for (k=0; k<l1; k++) {
            for (i=2; i<ido; i+=2) {
                  ic = ido - i;
                  VECLOOP {
                  tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido);
                  cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2;
                  ref(ch, i - 1 + k*ido) = ref(cc,i - 1 + 3*k*ido) + tr2;
                  ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido);
                  ci2 = ref(cc,i + 3*k*ido) + taur*ti2;
                  ref(ch, i + k*ido) = ref(cc,i + 3*k*ido) + ti2;
                  cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido));
                  ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido));
                  dr2 = cr2 - ci3;
                  dr3 = cr2 + ci3;
                  di2 = ci2 + cr3;
                  di3 = ci2 - cr3;
                  ref(ch, i - 1 + (k + l1)*ido) = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
                  ref(ch, i + (k + l1)*ido) = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
                  ref(ch, i - 1 + (k + 2*l1)*ido) = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
                  ref(ch, i + (k + 2*l1)*ido) = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
                  }
            }
    }
} /* radb3 */


static void radf4(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[], const Treal wa3[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
    for (k=0; k<l1; k++) {
            VECLOOP {
            tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido);
            tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido);
            ref(ch, 4*k*ido) = tr1 + tr2;
            ref(ch, ido-1 + (4*k + 3)*ido) = tr2 - tr1;
            ref(ch, ido-1 + (4*k + 1)*ido) = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido);
            ref(ch, (4*k + 2)*ido) = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido);
            }
    }
    if (ido < 2) return;
    if (ido != 2) {
            for (k=0; k<l1; k++) {
                  for (i=2; i<ido; i += 2) {
                        VECLOOP {
                        ic = ido - i;
                        cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
                        ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
                        cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
                        ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
                        cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
                        ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
                        tr1 = cr2 + cr4;
                        tr4 = cr4 - cr2;
                        ti1 = ci2 + ci4;
                        ti4 = ci2 - ci4;
                        ti2 = ref(cc,i + k*ido) + ci3;
                        ti3 = ref(cc,i + k*ido) - ci3;
                        tr2 = ref(cc,i - 1 + k*ido) + cr3;
                        tr3 = ref(cc,i - 1 + k*ido) - cr3;
                        ref(ch, i - 1 + 4*k*ido) = tr1 + tr2;
                        ref(ch, ic - 1 + (4*k + 3)*ido) = tr2 - tr1;
                        ref(ch, i + 4*k*ido) = ti1 + ti2;
                        ref(ch, ic + (4*k + 3)*ido) = ti1 - ti2;
                        ref(ch, i - 1 + (4*k + 2)*ido) = ti4 + tr3;
                        ref(ch, ic - 1 + (4*k + 1)*ido) = tr3 - ti4;
                        ref(ch, i + (4*k + 2)*ido) = tr4 + ti3;
                        ref(ch, ic + (4*k + 1)*ido) = tr4 - ti3;
                        }
                  }
            }
            if (ido % 2 == 1) return;
      }
    for (k=0; k<l1; k++) {
            VECLOOP {
            ti1 = -invsqrt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido));
            tr1 = invsqrt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido));
            ref(ch, ido-1 + 4*k*ido) = tr1 + ref(cc,ido-1 + k*ido);
            ref(ch, ido-1 + (4*k + 2)*ido) = ref(cc,ido-1 + k*ido) - tr1;
            ref(ch, (4*k + 1)*ido) = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido);
            ref(ch, (4*k + 3)*ido) = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido);
            }
    }
} /* radf4 */


static void radb4(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[], const Treal wa3[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
    for (k = 0; k < l1; k++) {
            VECLOOP {
            tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido);
            tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido);
            tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido);
            tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido);
            ref(ch, k*ido) = tr2 + tr3;
            ref(ch, (k + l1)*ido) = tr1 - tr4;
            ref(ch, (k + 2*l1)*ido) = tr2 - tr3;
            ref(ch, (k + 3*l1)*ido) = tr1 + tr4;
            }
    }
    if (ido < 2) return;
    if (ido != 2) {
            for (k = 0; k < l1; ++k) {
                  for (i = 2; i < ido; i += 2) {
                        ic = ido - i;
                        VECLOOP {
                        ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido);
                        ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido);
                        ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido);
                        tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido);
                        tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido);
                        tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido);
                        ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido);
                        tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido);
                        ref(ch, i - 1 + k*ido) = tr2 + tr3;
                        cr3 = tr2 - tr3;
                        ref(ch, i + k*ido) = ti2 + ti3;
                        ci3 = ti2 - ti3;
                        cr2 = tr1 - tr4;
                        cr4 = tr1 + tr4;
                        ci2 = ti1 + ti4;
                        ci4 = ti1 - ti4;
                        ref(ch, i - 1 + (k + l1)*ido) = wa1[i - 2]*cr2 - wa1[i - 1]*ci2;
                        ref(ch, i + (k + l1)*ido) = wa1[i - 2]*ci2 + wa1[i - 1]*cr2;
                        ref(ch, i - 1 + (k + 2*l1)*ido) = wa2[i - 2]*cr3 - wa2[i - 1]*ci3;
                        ref(ch, i + (k + 2*l1)*ido) = wa2[i - 2]*ci3 + wa2[i - 1]*cr3;
                        ref(ch, i - 1 + (k + 3*l1)*ido) = wa3[i - 2]*cr4 - wa3[i - 1]*ci4;
                        ref(ch, i + (k + 3*l1)*ido) = wa3[i - 2]*ci4 + wa3[i - 1]*cr4;
                        }
                  }
            }
            if (ido % 2 == 1) return;
      }
    for (k = 0; k < l1; k++) {
            VECLOOP {
            ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido);
            ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido);
            tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido);
            tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido);
            ref(ch, ido-1 + k*ido) = tr2 + tr2;
            ref(ch, ido-1 + (k + l1)*ido) = sqrt2*(tr1 - ti1);
            ref(ch, ido-1 + (k + 2*l1)*ido) = ti2 + ti2;
            ref(ch, ido-1 + (k + 3*l1)*ido) = -sqrt2*(tr1 + ti1);
            }
    }
} /* radb4 */


static void radf5(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5,
           cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
    for (k = 0; k < l1; k++) {
            VECLOOP {
            cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido);
            ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido);
            cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido);
            ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido);
            ref(ch, 5*k*ido) = ref(cc,k*ido) + cr2 + cr3;
            ref(ch, ido-1 + (5*k + 1)*ido) = ref(cc,k*ido) + tr11*cr2 + tr12*cr3;
            ref(ch, (5*k + 2)*ido) = ti11*ci5 + ti12*ci4;
            ref(ch, ido-1 + (5*k + 3)*ido) = ref(cc,k*ido) + tr12*cr2 + tr11*cr3;
            ref(ch, (5*k + 4)*ido) = ti12*ci5 - ti11*ci4;
            }
    }
    if (ido == 1) return;
    for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
                  VECLOOP {
                  ic = ido - i;
                  dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido);
                  di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido);
                  dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido);
                  di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido);
                  dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido);
                  di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido);
                  dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido);
                  di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido);
                  cr2 = dr2 + dr5;
                  ci5 = dr5 - dr2;
                  cr5 = di2 - di5;
                  ci2 = di2 + di5;
                  cr3 = dr3 + dr4;
                  ci4 = dr4 - dr3;
                  cr4 = di3 - di4;
                  ci3 = di3 + di4;
                  ref(ch, i - 1 + 5*k*ido) = ref(cc,i - 1 + k*ido) + cr2 + cr3;
                  ref(ch, i + 5*k*ido) = ref(cc,i + k*ido) + ci2 + ci3;
                  tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3;
                  ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3;
                  tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3;
                  ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3;
                  tr5 = ti11*cr5 + ti12*cr4;
                  ti5 = ti11*ci5 + ti12*ci4;
                  tr4 = ti12*cr5 - ti11*cr4;
                  ti4 = ti12*ci5 - ti11*ci4;
                  ref(ch, i - 1 + (5*k + 2)*ido) = tr2 + tr5;
                  ref(ch, ic - 1 + (5*k + 1)*ido) = tr2 - tr5;
                  ref(ch, i + (5*k + 2)*ido) = ti2 + ti5;
                  ref(ch, ic + (5*k + 1)*ido) = ti5 - ti2;
                  ref(ch, i - 1 + (5*k + 4)*ido) = tr3 + tr4;
                  ref(ch, ic - 1 + (5*k + 3)*ido) = tr3 - tr4;
                  ref(ch, i + (5*k + 4)*ido) = ti3 + ti4;
                  ref(ch, ic + (5*k + 3)*ido) = ti4 - ti3;
                  }
            }
    }
} /* radf5 */


static void radb5(Tint ido, Tint l1, const Treal cc[], Treal ch[], 
                          const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k, ic;
    Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
           ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
    for (k = 0; k < l1; k++) {
            VECLOOP {
            ti5 = 2*ref(cc,(5*k + 2)*ido);
            ti4 = 2*ref(cc,(5*k + 4)*ido);
            tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido);
            tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido);
            ref(ch, k*ido) = ref(cc,5*k*ido) + tr2 + tr3;
            cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3;
            cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3;
            ci5 = ti11*ti5 + ti12*ti4;
            ci4 = ti12*ti5 - ti11*ti4;
            ref(ch, (k + l1)*ido) = cr2 - ci5;
            ref(ch, (k + 2*l1)*ido) = cr3 - ci4;
            ref(ch, (k + 3*l1)*ido) = cr3 + ci4;
            ref(ch, (k + 4*l1)*ido) = cr2 + ci5;
            }
    }
    if (ido == 1) return;
    for (k = 0; k < l1; ++k) {
            for (i = 2; i < ido; i += 2) {
                  VECLOOP {
                  ic = ido - i;
                  ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido);
                  ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido);
                  ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido);
                  ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido);
                  tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido);
                  tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido);
                  tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido);
                  tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido);
                  ref(ch, i - 1 + k*ido) = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3;
                  ref(ch, i + k*ido) = ref(cc,i + 5*k*ido) + ti2 + ti3;
                  cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3;

                  ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3;
                  cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3;

                  ci3 = ref(cc,i + 5*k*ido) + 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;
                  ref(ch, i - 1 + (k + l1)*ido) = wa1[i - 2]*dr2 - wa1[i - 1]*di2;
                  ref(ch, i + (k + l1)*ido) = wa1[i - 2]*di2 + wa1[i - 1]*dr2;
                  ref(ch, i - 1 + (k + 2*l1)*ido) = wa2[i - 2]*dr3 - wa2[i - 1]*di3;
                  ref(ch, i + (k + 2*l1)*ido) = wa2[i - 2]*di3 + wa2[i - 1]*dr3;
                  ref(ch, i - 1 + (k + 3*l1)*ido) = wa3[i - 2]*dr4 - wa3[i - 1]*di4;
                  ref(ch, i + (k + 3*l1)*ido) = wa3[i - 2]*di4 + wa3[i - 1]*dr4;
                  ref(ch, i - 1 + (k + 4*l1)*ido) = wa4[i - 2]*dr5 - wa4[i - 1]*di5;
                  ref(ch, i + (k + 4*l1)*ido) = wa4[i - 2]*di5 + wa4[i - 1]*dr5;
                  }
            }
    }
} /* radb5 */


static void radfg(Tint ido, Tint ip, Tint l1, Tint idl1,
                          Treal cc[], Treal ch[], const Treal wa[] VECARGS_FORMAL)
{
      VECVARS;
    Tint idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
    Treal dcc, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
    arg = twopi / ip;
    dcp = cos(arg);
    dsp = sin(arg);
    ipph = (ip + 1) / 2;
    nbd = (ido - 1) / 2;
    if (ido != 1) {
            for (ik=0; ik<idl1; ik++) {VECLOOP ref(ch, ik) = ref(cc, ik);}
            for (j=1; j<ip; j++)
                  for (k=0; k<l1; k++) {
                        VECLOOP
                              ref(ch, (k + j*l1)*ido) = ref(cc, (k + j*l1)*ido);
                  }
            if (nbd <= l1) {
                  is = -ido;
                  for (j=1; j<ip; j++) {
                        is += ido;
                        idij = is-1;
                        for (i=2; i<ido; i+=2) {
                              idij += 2;
                              for (k=0; k<l1; k++) {
                                    VECLOOP {
                                    ref(ch, i - 1 + (k + j*l1)*ido) =
                                          wa[idij - 1]*ref(cc, i - 1 + (k + j*l1)*ido) + wa[idij]*ref(cc, i + (k + j*l1)*ido);
                                    ref(ch, i + (k + j*l1)*ido) =
                                          wa[idij - 1]*ref(cc, i + (k + j*l1)*ido) - wa[idij]*ref(cc, i - 1 + (k + j*l1)*ido);
                                    }
                              }
                        }
                  }
            } else {
                  is = -ido;
                  for (j=1; j<ip; j++) {
                        is += ido;
                        for (k=0; k<l1; k++) {
                              idij = is-1;
                              for (i=2; i<ido; i+=2) {
                                    idij += 2;
                                    VECLOOP {
                                    ref(ch, i - 1 + (k + j*l1)*ido) =
                                          wa[idij - 1]*ref(cc, i - 1 + (k + j*l1)*ido) + wa[idij]*ref(cc, i + (k + j*l1)*ido);
                                    ref(ch, i + (k + j*l1)*ido) =
                                          wa[idij - 1]*ref(cc, i + (k + j*l1)*ido) - wa[idij]*ref(cc, i - 1 + (k + j*l1)*ido);
                                    }
                              }
                        }
                  }
            }
            if (nbd >= l1) {
                  for (j=1; j<ipph; j++) {
                        jc = ip - j;
                        for (k=0; k<l1; k++) {
                              for (i=2; i<ido; i+=2) {
                                    VECLOOP {
                                    ref(cc, i - 1 + (k + j*l1)*ido) =
                                          ref(ch, i - 1 + (k + j*l1)*ido) + ref(ch, i - 1 + (k + jc*l1)*ido);
                                    ref(cc, i - 1 + (k + jc*l1)*ido) =
                                          ref(ch, i + (k + j*l1)*ido) - ref(ch, i + (k + jc*l1)*ido);
                                    ref(cc, i + (k + j*l1)*ido) =
                                          ref(ch, i + (k + j*l1)*ido) + ref(ch, i + (k + jc*l1)*ido);
                                    ref(cc, i + (k + jc*l1)*ido) =
                                          ref(ch, i - 1 + (k + jc*l1)*ido) - ref(ch, i - 1 + (k + j*l1)*ido);
                                    }
                              }
                        }
                  }
            } else {
                  for (j=1; j<ipph; j++) {
                        jc = ip - j;
                        for (i=2; i<ido; i+=2) {
                              for (k=0; k<l1; k++) {
                                    VECLOOP {
                                    ref(cc, i - 1 + (k + j*l1)*ido) =
                                          ref(ch, i - 1 + (k + j*l1)*ido) + ref(ch, i - 1 + (k + jc*l1)*ido);
                                    ref(cc, i - 1 + (k + jc*l1)*ido) =
                                          ref(ch, i + (k + j*l1)*ido) - ref(ch, i + (k + jc*l1)*ido);
                                    ref(cc, i + (k + j*l1)*ido) =
                                          ref(ch, i + (k + j*l1)*ido) + ref(ch, i + (k + jc*l1)*ido);
                                    ref(cc, i + (k + jc*l1)*ido) =
                                          ref(ch, i - 1 + (k + jc*l1)*ido) - ref(ch, i - 1 + (k + j*l1)*ido);
                                    }
                              }
                        }
                  }
            }
      } else {    /* now ido == 1 */
            for (ik=0; ik<idl1; ik++) {VECLOOP ref(cc, ik) = ref(ch, ik);}
      }
      for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ref(cc, (k + j*l1)*ido) = ref(ch, (k + j*l1)*ido) + ref(ch, (k + jc*l1)*ido);
                  ref(cc, (k + jc*l1)*ido) = ref(ch, (k + jc*l1)*ido) - ref(ch, (k + j*l1)*ido);
                  }
            }
    }

    ar1 = 1;
    ai1 = 0;
    for (l=1; l<ipph; l++) {
            lc = ip - l;
            ar1h = dcp*ar1 - dsp*ai1;
            ai1 = dcp*ai1 + dsp*ar1;
            ar1 = ar1h;
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP {
                  ref(ch, ik + l*idl1) = ref(cc, ik) + ar1*ref(cc, ik + idl1);
                  ref(ch, ik + lc*idl1) = ai1*ref(cc, ik + (ip-1)*idl1);
                  }
            }
            dcc = ar1;
            ds2 = ai1;
            ar2 = ar1;
            ai2 = ai1;
            for (j=2; j<ipph; j++) {
                  jc = ip - j;
                  ar2h = dcc*ar2 - ds2*ai2;
                  ai2 = dcc*ai2 + ds2*ar2;
                  ar2 = ar2h;
                  for (ik=0; ik<idl1; ik++) {
                        VECLOOP {
                        ref(ch, ik + l*idl1) += ar2*ref(cc, ik + j*idl1);
                        ref(ch, ik + lc*idl1) += ai2*ref(cc, ik + jc*idl1);
                        }
                  }
            }
    }
    for (j=1; j<ipph; j++)
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP
                        ref(ch, ik) += ref(cc, ik + j*idl1);
            }

      if (ido >= l1) {
            for (k=0; k<l1; k++) {
                  for (i=0; i<ido; i++) {
                        VECLOOP {
                        ref(cc,i + k*ip*ido) = ref(ch, i + k*ido);
                        }
                  }
            }
      } else {
            for (i=0; i<ido; i++) {
                  for (k=0; k<l1; k++) {
                        VECLOOP {
                        ref(cc,i + k*ip*ido) = ref(ch, i + k*ido);
                        }
                  }
            }
      }
    for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) =
                        ref(ch, (k + j*l1)*ido);
                  ref(cc,(j2 + k*ip)*ido) =
                        ref(ch, (k + jc*l1)*ido);
                  }
            }
    }
    if (ido == 1) return;
    if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  j2 = 2*j;
                  for (k=0; k<l1; k++) {
                        for (i=2; i<ido; i+=2) {
                              ic = ido - i;
                              VECLOOP {
                              ref(cc,i - 1 + (j2 + k*ip)*ido) =
                                    ref(ch, i - 1 + (k + j*l1)*ido) + ref(ch, i - 1 + (k + jc*l1)*ido);
                              ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) =
                                    ref(ch, i - 1 + (k + j*l1)*ido) - ref(ch, i - 1 + (k + jc*l1)*ido);
                              ref(cc,i + (j2 + k*ip)*ido) =
                                    ref(ch, i + (k + j*l1)*ido) + ref(ch, i + (k + jc*l1)*ido);
                              ref(cc,ic + (j2 - 1 + k*ip)*ido) =
                                    ref(ch, i + (k + jc*l1)*ido) - ref(ch, i + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      } else {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  j2 = 2*j;
                  for (i=2; i<ido; i+=2) {
                        ic = ido - i;
                        for (k=0; k<l1; k++) {
                              VECLOOP {
                              ref(cc,i - 1 + (j2 + k*ip)*ido) = ref(ch, i - 1 + (k + j*l1)*ido) + ref(ch, i - 1 + (k + jc*l1)*ido);
                              ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ref(ch, i - 1 + (k + j*l1)*ido) - ref(ch, i - 1 + (k + jc*l1)*ido);
                              ref(cc,i + (j2 + k*ip)*ido) = ref(ch, i + (k + j*l1)*ido) + ref(ch, i + (k + jc*l1)*ido);
                              ref(cc,ic + (j2 - 1 + k*ip)*ido) = ref(ch, i + (k + jc*l1)*ido) - ref(ch, i + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      }
} /* radfg */


static void radbg(Tint ido, Tint ip, Tint l1, Tint idl1,
                          Treal cc[], Treal ch[], const Treal wa[] VECARGS_FORMAL)
{
      VECVARS;
    Tint idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
    Treal dcc, ai1, ai2, ar1, ar2, ds2;
    Tint nbd;
    Treal dcp, arg, dsp, ar1h, ar2h;
    arg = twopi / ip;
    dcp = cos(arg);
    dsp = sin(arg);
    nbd = (ido - 1) / 2;
    ipph = (ip + 1) / 2;
    if (ido >= l1) {
            for (k=0; k<l1; k++) {
                  for (i=0; i<ido; i++) {
                        VECLOOP {
                        ref(ch, i + k*ido) = ref(cc,i + k*ip*ido);
                        }
                  }
            }
      } else {
            for (i=0; i<ido; i++) {
                  for (k=0; k<l1; k++) {
                        VECLOOP {
                        ref(ch, i + k*ido) = ref(cc,i + k*ip*ido);
                        }
                  }
            }
      }
    for (j=1; j<ipph; j++) {
            jc = ip - j;
            j2 = 2*j;
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ref(ch, (k + j*l1)*ido) = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)*ido);
                  ref(ch, (k + jc*l1)*ido) = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido);
                  }
            }
    }

    if (ido != 1) {
            if (nbd >= l1) {
                  for (j=1; j<ipph; j++) {
                        jc = ip - j;
                        for (k=0; k<l1; k++) {
                              for (i=2; i<ido; i+=2) {
                                    ic = ido - i;
                                    VECLOOP {
                                    ref(ch, i - 1 + (k + j*l1)*ido) =
                                          ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i - 1 + (k + jc*l1)*ido) =
                                          ref(cc,i - 1 + (2*j + k*ip)*ido) - ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i + (k + j*l1)*ido) =
                                          ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i + (k + jc*l1)*ido) =
                                          ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic + (2*j - 1 + k*ip)*ido);
                                    }
                              }
                        }
                  }
            } else {
                  for (j=1; j<ipph; j++) {
                        jc = ip - j;
                        for (i=2; i<ido; i+=2) {
                              ic = ido - i;
                              for (k=0; k<l1; k++) {
                                    VECLOOP {
                                    ref(ch, i - 1 + (k + j*l1)*ido) =
                                          ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i - 1 + (k + jc*l1)*ido) =
                                          ref(cc,i - 1 + (2*j + k*ip)*ido) - ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i + (k + j*l1)*ido) =
                                          ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic + (2*j - 1 + k*ip)*ido);
                                    ref(ch, i + (k + jc*l1)*ido) =
                                          ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic + (2*j - 1 + k*ip)*ido);
                                    }
                              }
                        }
                  }
            }
      }

      ar1 = 1;
    ai1 = 0;
    for (l=1; l<ipph; l++) {
            lc = ip - l;
            ar1h = dcp*ar1 - dsp*ai1;
            ai1 = dcp*ai1 + dsp*ar1;
            ar1 = ar1h;
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP {
                  ref(cc, ik + l*idl1) = ref(ch, ik) + ar1*ref(ch, ik + idl1);
                  ref(cc, ik + lc*idl1) = ai1*ref(ch, ik + (ip-1)*idl1);
                  }
            }
            dcc = ar1;
            ds2 = ai1;
            ar2 = ar1;
            ai2 = ai1;
            for (j=2; j<ipph; j++) {
                  jc = ip - j;
                  ar2h = dcc*ar2 - ds2*ai2;
                  ai2 = dcc*ai2 + ds2*ar2;
                  ar2 = ar2h;
                  for (ik=0; ik<idl1; ik++) {
                        VECLOOP {
                        ref(cc, ik + l*idl1) += ar2*ref(ch, ik + j*idl1);
                        ref(cc, ik + lc*idl1) += ai2*ref(ch, ik + jc*idl1);
                        }
                  }
            }
    }
    for (j=1; j<ipph; j++) {
            for (ik=0; ik<idl1; ik++) {
                  VECLOOP {
                  ref(ch, ik) += ref(ch, ik + j*idl1);
                  }
            }
    }
    for (j=1; j<ipph; j++) {
            jc = ip - j;
            for (k=0; k<l1; k++) {
                  VECLOOP {
                  ref(ch, (k + j*l1)*ido) = ref(cc, (k + j*l1)*ido) - ref(cc, (k + jc*l1)*ido);
                  ref(ch, (k + jc*l1)*ido) = ref(cc, (k + j*l1)*ido) + ref(cc, (k + jc*l1)*ido);
                  }
            }
    }
      
    if (ido == 1) return;
      if (nbd >= l1) {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  for (k=0; k<l1; k++) {
                        for (i=2; i<ido; i+=2) {
                              VECLOOP {
                              ref(ch, i - 1 + (k + j*l1)*ido) = ref(cc, i - 1 + (k + j*l1)*ido) - ref(cc, i + (k + jc*l1)*ido);
                              ref(ch, i - 1 + (k + jc*l1)*ido) = ref(cc, i - 1 + (k + j*l1)*ido) + ref(cc, i + (k + jc*l1)*ido);
                              ref(ch, i + (k + j*l1)*ido) = ref(cc, i + (k + j*l1)*ido) + ref(cc, i - 1 + (k + jc*l1)*ido);
                              ref(ch, i + (k + jc*l1)*ido) = ref(cc, i + (k + j*l1)*ido) - ref(cc, i - 1 + (k + jc*l1)*ido);
                              }
                        }
                  }
            }
      } else {
            for (j=1; j<ipph; j++) {
                  jc = ip - j;
                  for (i=2; i<ido; i+=2) {
                        for (k=0; k<l1; k++) {
                              VECLOOP {
                              ref(ch, i - 1 + (k + j*l1)*ido) = ref(cc, i - 1 + (k + j*l1)*ido) - ref(cc, i + (k + jc*l1)*ido);
                              ref(ch, i - 1 + (k + jc*l1)*ido) = ref(cc, i - 1 + (k + j *l1)*ido) + ref(cc, i + (k + jc*l1)*ido);
                              ref(ch, i + (k + j*l1)*ido) = ref(cc, i + (k + j*l1)*ido) + ref(cc, i - 1 + (k + jc*l1)*ido);
                              ref(ch, i + (k + jc*l1)*ido) = ref(cc, i + (k + j*l1)*ido) - ref(cc, i - 1 + (k + jc*l1)*ido);
                              }
                        }
                  }
            }
      }
    for (ik=0; ik<idl1; ik++) {VECLOOP ref(cc, ik) = ref(ch, ik);}
    for (j=1; j<ip; j++)
            for (k=0; k<l1; k++) {
                  VECLOOP
                        ref(cc, (k + j*l1)*ido) = ref(ch, (k + j*l1)*ido);
            }
    if (nbd <= l1) {
            is = -ido;
            for (j=1; j<ip; j++) {
                  is += ido;
                  idij = is-1;
                  for (i=2; i<ido; i+=2) {
                        idij += 2;
                        for (k=0; k<l1; k++) {
                              VECLOOP {
                              ref(cc, i - 1 + (k + j*l1)*ido) =
                                    wa[idij - 1]*ref(ch, i - 1 + (k + j*l1)*ido) - wa[idij]*ref(ch, i + (k + j*l1)*ido);
                              ref(cc, i + (k + j*l1)*ido) =
                                    wa[idij - 1]*ref(ch, i + (k + j*l1)*ido) + wa[idij]*ref(ch, i - 1 + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      } else {
            is = -ido;
            for (j=1; j<ip; j++) {
                  is += ido;
                  for (k=0; k<l1; k++) {
                        idij = is;
                        for (i=2; i<ido; i+=2) {
                              idij += 2;
                              VECLOOP {
                              ref(cc, i - 1 + (k + j*l1)*ido) =
                                    wa[idij-1]*ref(ch, i - 1 + (k + j*l1)*ido) - wa[idij]*ref(ch, i + (k + j*l1)*ido);
                              ref(cc, i + (k + j*l1)*ido) =
                                    wa[idij-1]*ref(ch, i + (k + j*l1)*ido) + wa[idij]*ref(ch, i - 1 + (k + j*l1)*ido);
                              }
                        }
                  }
            }
      }
} /* radbg */


/* ----------------------------------------------------------------------
   cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
   ---------------------------------------------------------------------- */

static void cfftf1(Tint n, Treal c[], Treal ch[], const Treal wa[], const Tint ifac[MAXFAC+2], int isign VECARGS_FORMAL)
{
      VECVARS;
    Tint idot, i;
    Tint k1, l1=1, l2;
    Tint na=0, nf=ifac[1], ip, iw=0, ix2, ix3, ix4, nac, ido, idl1;
      Treal *cinput, *coutput;
    for (k1=2; k1<=nf+1; k1++) {
            ip = ifac[k1];
            l2 = ip*l1;
            ido = n / l2;
            idot = 2*ido;
            idl1 = idot*l1;
            if (na) {
                  cinput = ch;
                  coutput = c;
            } else {
                  cinput = c;
                  coutput = ch;
            }
            switch (ip) {
             case 4:
                  ix2 = iw + idot;
                  ix3 = ix2 + idot;
                  passf4(idot, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3, isign VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 2:
                  passf2(idot, l1, cinput, coutput, wa+iw, isign VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 3:
                  ix2 = iw + idot;
                  passf3(idot, l1, cinput, coutput, wa+iw, wa+ix2, isign VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 5:
                  ix2 = iw + idot;
                  ix3 = ix2 + idot;
                  ix4 = ix3 + idot;
                  passf5(idot, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3, wa+ix4, isign VECARGS_ACTUAL);
                  na = !na;
                  break;
             default:
                  passf(&nac, idot, ip, l1, idl1, cinput, coutput, wa+iw, isign VECARGS_ACTUAL);
                  if (nac) na = !na;
            }
            l1 = l2;
            iw += (ip - 1)*idot;
    }
      if (na)
            for (i=0; i<2*n; i++) {VECLOOP ref0(c,i) = ref0(ch,i);}
} /* cfftf1 */


void veccfftf(Tint n, Treal c[], Treal work[] VECARGS_FORMAL)
/* Complex forward transform.
   Work array length at least 2*((lot-1)*jump + n*step) + 2*n + MAXFAC + 2. */
{
    Tint iw1, iw2;
    if (n <= 1) return;
    iw1 = 2*((lot-1)*jump + n*step);
    iw2 = iw1 + 2*n;
    cfftf1(n, c, work, work+iw1, (Tint*)(work+iw2), -1 VECARGS_ACTUAL);
} /* cfftf */


void veccfftb(Tint n, Treal c[], Treal work[] VECARGS_FORMAL)
/* Complex backward transform.
   Work array length at least 2*((lot-1)*jump + n*step) + 2*n + MAXFAC + 2. */
{
    Tint iw1, iw2;
    if (n <= 1) return;
    iw1 = 2*((lot-1)*jump + n*step);
    iw2 = iw1 + 2*n;
    cfftf1(n, c, work, work+iw1, (Tint*)(work+iw2), +1 VECARGS_ACTUAL);
} /* cfftb */


static void factorize(Tint n, Tint ifac[MAXFAC+2])
/* Factorize n in factors of 2,3,4,5 and rest. On exit,
   ifac[0] contains n and ifac[1] contains number of factors,
   the factors start from ifac[2]. */
{
    static const Tint ntryh[4] = {3,4,2,5};
    Tint ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr;
      if (sizeof(Tint) > sizeof(Treal)) {
            fprintf(stderr,"*** vecfftpack: sizeof(Tint) > sizeof(Treal)\n");
            exit(1);
      }
 startloop:
      if (j < 4)
            ntry = ntryh[j];
      else
            ntry+= 2;
      j++;
      do {
            nq = nl / ntry;
            nr = nl - ntry*nq;
            if (nr != 0) goto startloop;
            nf++;
            ifac[nf + 1] = ntry;
            nl = nq;
            if (ntry == 2 && nf != 1) {
                  for (i=2; i<=nf; i++) {
                        ib = nf - i + 2;
                        ifac[ib + 1] = ifac[ib];
                  }
                  ifac[2] = 2;
            }
      } while (nl != 1);
    ifac[0] = n;
    ifac[1] = nf;
}


static void cffti1(Tint n, Treal wa[], Tint ifac[MAXFAC+2])
{
    Treal arg, argh, argld, fi;
    Tint ido, i, j, idot, ipm;
    Tint i1, k1, l1, l2;
    Tint ld, ii, nf, ip;
      factorize(n,ifac);
      nf = ifac[1];
    argh = twopi/(Treal)n;
    i = 1;
    l1 = 1;
    for (k1=1; k1<=nf; k1++) {
            ip = ifac[k1+1];
            ld = 0;
            l2 = l1*ip;
            ido = n / l2;
            idot = ido + ido + 2;
            ipm = ip - 1;
            for (j=1; j<=ipm; j++) {
                  i1 = i;
                  wa[i-1] = 1;
                  wa[i] = 0;
                  ld += l1;
                  fi = 0;
                  argld = ld*argh;
                  for (ii=4; ii<=idot; ii+=2) {
                        i+= 2;
                        fi+= 1;
                        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;
    }
} /* cffti1 */


void veccffti(Tint n, Treal work[] VECARGS_FORMAL)
/* Initialize complex forward or backward transform.
   Work array length at least 2*((lot-1)*jump + n*step) + 2*n + MAXFAC + 2. */
{
    Tint iw1, iw2;
    if (n <= 1) return;
      iw1 = 2*((lot-1)*jump + n*step);
    iw2 = iw1 + 2*n;
    cffti1(n, work+iw1, (Tint*)(work+iw2));
} /* cffti */

/* ----------------------------------------------------------------------
   rfftf1, rfftb1, rfftf, rfftb, rffti1, rffti. Real FFTs.
   ---------------------------------------------------------------------- */

static void rfftf1(Tint n, Treal c[], Treal ch[], const Treal wa[], const Tint ifac[MAXFAC+2] VECARGS_FORMAL)
{
      VECVARS;
    Tint i;
    Tint k1, l1, l2=n, na=1, kh, nf=ifac[1], ip, iw=n-1, ix2, ix3, ix4, ido, idl1;
      Treal *cinput, *coutput;
    for (k1=1; k1<=nf; k1++) {
            kh = nf - k1;
            ip = ifac[kh + 2];
            l1 = l2 / ip;
            ido = n / l2;
            idl1 = ido*l1;
            iw -= (ip - 1)*ido;
            na = !na;
            if (na) {
                  cinput = ch;
                  coutput = c;
            } else {
                  cinput = c;
                  coutput = ch;
            }
            switch (ip) {
             case 4:
                  ix2 = iw + ido;
                  ix3 = ix2 + ido;
                  radf4(ido, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3 VECARGS_ACTUAL);
                  break;
             case 2:
                  radf2(ido, l1, cinput, coutput, wa+iw VECARGS_ACTUAL);
                  break;
             case 3:
                  ix2 = iw + ido;
                  radf3(ido, l1, cinput, coutput, wa+iw, wa+ix2 VECARGS_ACTUAL);
                  break;
             case 5:
                  ix2 = iw + ido;
                  ix3 = ix2 + ido;
                  ix4 = ix3 + ido;
                  radf5(ido, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3, wa+ix4 VECARGS_ACTUAL);
                  break;
             default:
                  if (ido == 1) na = !na;
                  if (na)
                        radfg(ido, ip, l1, idl1, ch, c, wa+iw VECARGS_ACTUAL);
                  else
                        radfg(ido, ip, l1, idl1, c, ch, wa+iw VECARGS_ACTUAL);
                  na = !na;
            }
            l2 = l1;
    }
    if (na) return;
    for (i=0; i<n; i++) {VECLOOP ref(c,i) = ref(ch,i);}
} /* rfftf1 */


static void rfftb1(Tint n, Treal c[], Treal ch[], const Treal wa[], const Tint ifac[MAXFAC+2] VECARGS_FORMAL)
{
      VECVARS;
    Tint i, k1, l1=1, l2, na=0, nf=ifac[1], ip, iw=0, ix2, ix3, ix4, ido, idl1;
      Treal *cinput, *coutput;
    for (k1=1; k1<=nf; k1++) {
            ip = ifac[k1 + 1];
            l2 = ip*l1;
            ido = n / l2;
            idl1 = ido*l1;
            if (na) {
                  cinput = ch;
                  coutput = c;
            } else {
                  cinput = c;
                  coutput = ch;
            }
            switch (ip) {
             case 4:
                  ix2 = iw + ido;
                  ix3 = ix2 + ido;
                  radb4(ido, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3 VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 2:
                  radb2(ido, l1, cinput, coutput, wa+iw VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 3:
                  ix2 = iw + ido;
                  radb3(ido, l1, cinput, coutput, wa+iw, wa+ix2 VECARGS_ACTUAL);
                  na = !na;
                  break;
             case 5:
                  ix2 = iw + ido;
                  ix3 = ix2 + ido;
                  ix4 = ix3 + ido;
                  radb5(ido, l1, cinput, coutput, wa+iw, wa+ix2, wa+ix3, wa+ix4 VECARGS_ACTUAL);
                  na = !na;
                  break;
             default:
                  radbg(ido, ip, l1, idl1, cinput, coutput, wa+iw VECARGS_ACTUAL);
                  if (ido == 1) na = !na;
            }
            l1 = l2;
            iw += (ip - 1)*ido;
    }
    if (na)
            for (i=0; i<n; i++) {VECLOOP ref(c,i) = ref(ch,i);}
} /* rfftb1 */


void vecrfftf(Tint n, Treal r[], Treal work[] VECARGS_FORMAL)
/* Real forward transform.
   Work array length at least (lot-1)*jump + n*step + n + MAXFAC + 2. */
{
    if (n <= 1) return;
    rfftf1(n, r, work, work+(lot-1)*jump+n*step, (Tint*)(work+(lot-1)*jump+n*step+n) VECARGS_ACTUAL);
} /* rfftf */


void vecrfftb(Tint n, Treal r[], Treal work[] VECARGS_FORMAL)
/* Real backward transform.
   Work array length at least (lot-1)*jump + n*step + n + MAXFAC + 2. */
{
    if (n <= 1) return;
    rfftb1(n, r, work, work+(lot-1)*jump+n*step, (Tint*)(work+(lot-1)*jump+n*step+n) VECARGS_ACTUAL);
} /* rfftb */


static void rffti1(Tint n, Treal wa[], Tint ifac[MAXFAC+2])
{
    Treal arg, argh, argld, fi;
    Tint i, j;
    Tint k1, l1, l2;
    Tint ld, ii, nf, ip, is;
    Tint ido, ipm;
    Tint nfm1;
      factorize(n,ifac);
      nf = ifac[1];
    argh = twopi / (Treal) (n);
    is = 0;
    nfm1 = nf - 1;
    l1 = 1;
    if (nfm1 == 0) return;
    for (k1 = 1; k1 <= nfm1; k1++) {
            ip = ifac[k1 + 1];
            ld = 0;
            l2 = l1*ip;
            ido = n / l2;
            ipm = ip - 1;
            for (j = 1; j <= ipm; ++j) {
                  ld += l1;
                  i = is;
                  argld = (Treal) ld*argh;
                  fi = 0;
                  for (ii = 3; ii <= ido; ii += 2) {
                        i += 2;
                        fi += 1;
                        arg = fi*argld;
                        wa[i - 2] = cos(arg);
                        wa[i - 1] = sin(arg);
                  }
                  is += ido;
            }
            l1 = l2;
    }
} /* rffti1 */


void vecrffti(Tint n, Treal work[] VECARGS_FORMAL)
/* Initialize real forward or backward transform.
   Work array length at least (lot-1)*jump + n*step + n + MAXFAC + 2. */
{
    if (n <= 1) return;
    rffti1(n, work+(lot-1)*jump+n*step, (Tint*)(work+(lot-1)*jump+n*step+n));
} /* rffti */

/* ----------------------------------------------------------------------
   cosqf1, cosqb1, cosqf, cosqb, cosqi. Quarter-wave cos-FFTs.
   ---------------------------------------------------------------------- */

static void cosqf1(Tint n, Treal x[], Treal w[], Treal xh[] VECARGS_FORMAL)
{
      VECVARS;
    Tint modn, i, k;
    Tint kc, np2, ns2;
    Treal xim1;
    ns2 = (n + 1) / 2;
    np2 = n + 2;
    for (k=1; k<ns2; k++) {
            kc = n - k;
            VECLOOP {
            ref(xh,k) = ref(x,k) + ref(x,kc);
            ref(xh,kc) = ref(x,k) - ref(x,kc);
            }
    }
    modn = n % 2;
    if (modn == 0) {
            VECLOOP ref(xh,ns2) = ref(x,ns2) + ref(x,ns2);
      }
    for (k=1; k<ns2; k++) {
            kc = n - k;
            VECLOOP {
            ref(x,k) = w[k - 1]*ref(xh,kc) + w[kc - 1]*ref(xh,k);
            ref(x,kc) = w[k - 1]*ref(xh,k) - w[kc - 1]*ref(xh,kc);
            }
    }
    if (modn == 0) {
            VECLOOP ref(x,ns2) = w[ns2-1]*ref(xh,ns2);
      }
    vecrfftf(n, x, xh VECARGS_ACTUAL);
    for (i=2; i<n; i+=2) {
            VECLOOP {
            xim1 = ref(x,i - 1) - ref(x,i);
            ref(x,i) = ref(x,i - 1) + ref(x,i);
            ref(x,i - 1) = xim1;
            }
    }
} /* cosqf1 */


static void cosqb1(Tint n, Treal x[], Treal w[], Treal xh[] VECARGS_FORMAL)
{
      VECVARS;
    Tint modn, i, k;
    Tint kc, ns2;
    Treal xim1;
    ns2 = (n + 1) / 2;
    for (i=2; i<n; i+=2) {
            VECLOOP {
            xim1 = ref(x,i - 1) + ref(x,i);
            ref(x,i) -= ref(x,i - 1);
            ref(x,i - 1) = xim1;
            }
    }
    VECLOOP ref(x,0) += ref(x,0);
    modn = n % 2;
    if (modn == 0) {
            VECLOOP ref(x,n-1) += ref(x,n-1);
      }
    vecrfftb(n, x, xh VECARGS_ACTUAL);
    for (k=1; k<ns2; k++) {
            kc = n - k;
            VECLOOP {
            ref(xh,k) = w[k - 1]*ref(x,kc) + w[kc - 1]*ref(x,k);
            ref(xh,kc) = w[k - 1]*ref(x,k) - w[kc - 1]*ref(x,kc);
            }
    }
    if (modn == 0) {
            VECLOOP ref(x,ns2) = w[ns2-1]*(ref(x,ns2) + ref(x,ns2));
      }
    for (k=1; k<ns2; k++) {
            kc = n - k;
            VECLOOP {
            ref(x,k) = ref(xh,k) + ref(xh,kc);
            ref(x,kc) = ref(xh,k) - ref(xh,kc);
            }
    }
      VECLOOP ref(x,0) += ref(x,0);
} /* cosqb1 */


void veccosqf(Tint n, Treal x[], Treal work[] VECARGS_FORMAL)
/* Quarter-wave forward cosine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
      VECVARS;
    Treal tsqx;
    if (n < 2) {
            return;
    } else if (n == 2) {
            VECLOOP {
            tsqx = sqrt2*ref(x,1);
            ref(x,1) = ref(x,0) - tsqx;
            ref(x,0) += tsqx;
          }
      }
    else
            cosqf1(n, x, work, work+(lot-1)*jump+n*step VECARGS_ACTUAL);
} /* cosqf */


void veccosqb(Tint n, Treal x[], Treal work[] VECARGS_FORMAL)
/* Quarter-wave backward cosine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
      VECVARS;
    Treal x1;
    if (n < 2) {
            VECLOOP ref(x,0) *= 4;
    } else if (n == 2) {
            VECLOOP {
            x1 = 4*(ref(x,0) + ref(x,1));
            ref(x,1) = tsqrt2*(ref(x,0) - ref(x,1));
            ref(x,0) = x1;
            }
    } else
            cosqb1(n, x, work, work+(lot-1)*jump+n*step VECARGS_ACTUAL);
} /* cosqb */


void veccosqi(Tint n, Treal work[] VECARGS_FORMAL)
/* Initialize quarter-wave forward or backward cosine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
    Tint k;
    Treal dt;
    dt = 0.5*pi/n;
    for (k=0; k<n; k++)
            work[k] = cos((k+1)*dt);
    vecrffti(n, work+(lot-1)*jump+n*step VECARGS_ACTUAL);
} /* cosqi */

/* ----------------------------------------------------------------------
   cost, costi. Full-wave cos-FFTs.
   ---------------------------------------------------------------------- */

void veccost(Tint n, Treal x[], Treal wsave[] VECARGS_FORMAL)
/* Forward/backward cosine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
      VECVARS;
/*    const int max_lot = 1024;*/
#define max_lot 2048
    Tint modn, i, k;
    Treal t1, t2;
    Tint kc;
    Treal xi;
    Tint nm1;
    Treal x1h;
    Tint ns2;
    Treal tx2, x1p3;
      Treal c1[max_lot], xim2[max_lot];
      if (lot > max_lot) {
            fprintf(stderr,"*** veccost: lot=%d > %d. Increase max_lot in vecfftpack.c:veccost().\n",lot,max_lot);
            return;
      }
    nm1 = n - 1;
    ns2 = n / 2;
    if (n - 2 < 0)
            return;
    else if (n == 2) {
            VECLOOP {
            x1h = ref(x,0) + ref(x,1);
            ref(x,1) = ref(x,0) - ref(x,1);
            ref(x,0) = x1h;
            }
    } else if (n == 3) {
            VECLOOP {
            x1p3 = ref(x,0) + ref(x,2);
            tx2 = ref(x,1) + ref(x,1);
            ref(x,1) = ref(x,0) - ref(x,2);
            ref(x,0) = x1p3 + tx2;
            ref(x,2) = x1p3 - tx2;
            }
      } else {
            VECLOOP {
            c1[p] = ref(x,0) - ref(x,n-1);
            ref(x,0) += ref(x,n-1);
            }
            for (k=1; k<ns2; k++) {
                  kc = nm1 - k;
                  VECLOOP {
                  t1 = ref(x,k) + ref(x,kc);
                  t2 = ref(x,k) - ref(x,kc);
                  c1[p] += wsave[kc]*t2;
                  t2 = wsave[k]*t2;
                  ref(x,k) = t1 - t2;
                  ref(x,kc) = t1 + t2;
                  }
            }
            modn = n % 2;
            if (modn != 0) {
                  VECLOOP ref(x,ns2) *= 2;
            }
            vecrfftf(nm1, x, wsave+(lot-1)*jump+n*step VECARGS_ACTUAL);
            VECLOOP {
            xim2[p] = ref(x,1);
            ref(x,1) = c1[p];
          }
            for (i=3; i<n; i+=2) {
                  VECLOOP {
                  xi = ref(x,i);
                  ref(x,i) = ref(x,i - 2) - ref(x,i - 1);
                  ref(x,i - 1) = xim2[p];
                  xim2[p] = xi;
                  }
            }
            if (modn != 0) {VECLOOP ref(x,n-1) = xim2[p];}
      }
} /* cost */


void veccosti(Tint n, Treal wsave[] VECARGS_FORMAL)
/* Initialize forward/backward cosine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
    Tint k,kc,ns2;
    Treal dt;
    if (n <= 3) return;
    ns2 = n / 2;
    dt = pi / (n-1);
    for (k = 1; k < ns2; k++) {
            kc = n - k - 1;
            wsave[k] = 2*sin(k*dt);
            wsave[kc] = 2*cos(k*dt);
    }
    vecrffti(n-1, wsave+(lot-1)*jump+n*step VECARGS_ACTUAL);
} /* costi */


/* ----------------------------------------------------------------------
   sinqf, sinqb, sinqi. Quarter-wave sin-FFTs. These call cosqf etc.
   ---------------------------------------------------------------------- */

void vecsinqf(Tint n, Treal x[], Treal wsave[] VECARGS_FORMAL)
/* Quarter-wave forward sine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
      VECVARS;
    Tint k;
    Treal xhold;
    Tint kc, ns2;
    if (n <= 1) return;
    ns2 = n / 2;
    for (k=0; k<ns2; k++) {
            kc = n - k - 1;
            VECLOOP {
            xhold = ref(x,k);
            ref(x,k) = ref(x,kc);
            ref(x,kc) = xhold;
            }
    }
    veccosqf(n, x, wsave VECARGS_ACTUAL);
    for (k=1; k<n; k+=2) {
            VECLOOP
                  ref(x,k) = -ref(x,k);
      }
} /* sinqf */


void vecsinqb(Tint n, Treal x[], Treal wsave[] VECARGS_FORMAL)
/* Quarter-wave backward sine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
      VECVARS;
    Tint k;
    Treal xhold;
    Tint kc, ns2;
    if (n <= 1) {
            VECLOOP ref(x,0) *= 4;
            return;
      }
    ns2 = n / 2;
    for (k=1; k<n; k+=2) {VECLOOP ref(x,k) = -ref(x,k);}
    veccosqb(n, x, wsave VECARGS_ACTUAL);
    for (k=0; k<ns2; k++) {
            kc = n - k - 1;
            VECLOOP {
            xhold = ref(x,k);
            ref(x,k) = ref(x,kc);
            ref(x,kc) = xhold;
            }
    }
} /* sinqb */


void vecsinqi(Tint n, Treal wsave[] VECARGS_FORMAL)
/* Initialize quarter-wave forward or backward sine transform.
   Work array length at least 2*((lot-1)*jump + n*step) + n + MAXFAC + 2. */
{
    veccosqi(n, wsave VECARGS_ACTUAL);
}

/* ----------------------------------------------------------------------
   sint1, sint, sinti. Full-wave sin-FFTs.
   ---------------------------------------------------------------------- */

#define refnp1(u,a) u[p*jump1+(a)*step]
            
static void sint1(Tint n, Treal war[], Treal war2[], const Treal was[], Treal xh[], const Treal x[],
                          const Tint ifac[MAXFAC+2] VECARGS_FORMAL)
/* Dimensions:
   was[n/2], x[n+1], xh and war2 are full grids with n+1, war is full grid with n */
{
      VECVARS;
    Tint modn, i, k;
    Treal xhold, t1, t2;
    Tint kc, ns2;
      const Tint jump1 = ((step==1) ? (jump+1) : (jump));
    if (n == 1) {
            VECLOOP ref(war,0) *= 2;
    } else if (n == 2) {
            VECLOOP {
            xhold = sqrt3*(ref(war,0) + ref(war,1));
            ref(war,1) = sqrt3*(ref(war,0) - ref(war,1));
            ref(war,0) = xhold;
            }
    } else {
            ns2 = n / 2;
            VECLOOP refnp1(xh,0) = 0;
            for (k=0; k<ns2; k++) {
                  kc = n - k - 1;
                  VECLOOP {
                  t1 = ref(war,k) - ref(war,kc);
                  t2 = was[k]*(ref(war,k) + ref(war,kc));
                  refnp1(xh,k + 1) = t1 + t2;
                  refnp1(xh,kc + 1) = t2 - t1;
                  }
            }
            modn = n % 2;
            if (modn != 0) {
                  VECLOOP
                        refnp1(xh,ns2 + 1) = 4*ref(war,ns2);
            }
            for (i=0; i<n; i++) {VECLOOP refnp1(war2,i) = ref(war,i);}
            rfftf1(n+1, xh, war2, x, ifac, lot,step,step==1?jump+1:jump);
            VECLOOP ref(war,0) = 0.5*refnp1(xh,0);
            for (i=2; i<n; i+=2) {
                  VECLOOP {
                  ref(war,i-1) = -refnp1(xh,i);
                  ref(war,i) = ref(war,i-2) + refnp1(xh,i-1);
                  }
            }
            if (modn == 0) {VECLOOP ref(war,n-1) = -refnp1(xh,n);}
      }
} /* sint1 */


void vecsint(Tint n, Treal x[], Treal work[] VECARGS_FORMAL)
/* Forward/backward sine transform.
   Work array length at least 2*((lot-1)*(jump+1) + (n+1)*step) + n + n/2 + 2 + MAXFAC + 2. */
{
    Tint iw1, iw2, iw3, iw4;
    iw1 = n/2;
    iw2 = iw1 + (lot-1)*(step==1?jump+1:jump) + (n + 1)*step;
      iw3 = iw2 + (lot-1)*(step==1?jump+1:jump) + (n + 1)*step;
    iw4 = iw3 + n + 1;
    sint1(n, x, work+iw1, work, work+iw2, work+iw3, (Tint*)(work+iw4) VECARGS_ACTUAL);
} /* sint */


void vecsinti(Tint n, Treal work[] VECARGS_FORMAL)
/* Initialize forward/backward sine transform.
   Work array length at least 2*((lot-1)*(jump+1) + (n+1)*step) + n + n/2 + 2 + MAXFAC + 2. */
{
    Tint k;
    Treal dt;
    if (n <= 1) return;
    dt = pi / (n+1);
    for (k=0; k<n/2; k++)
            work[k] = 2*sin((k+1)*dt);
    vecrffti(n+1, work+n/2 + (lot-1)*(step==1?jump+1:jump) + (n+1)*step, lot, step, step==1?jump+1:jump);
} /* sinti */

#ifdef __cplusplus
}
#endif

Generated by  Doxygen 1.6.0   Back to index