ComposerDesktopProject/dev/externals/fastconv/mxfftd.c

/*  
    mxfftd.c: double precision version of mxfft.c
	 
	This file is part of the CDP System.

    The CDP System is free software; you can redistribute it
    and/or modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    The CDP System is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with the CDP System; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA

*/

/* This program converted from the FORTRAN routines by Singleton in
 * Section 1.4 of  "Programs for Digital Signal Processing", IEEE Press, 1979.
 *  Conversion by Trevor Wishart and Keith Henderson, York Univ.
 */




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


void fft_(double *, double *, int, int, int, int);
void fftmx(double *, double *, int, int, int, int, int,
      int*, double *, double *, double *, double *, int *, int[]);
void reals_(double *, double *, int, int);

/*
 *-----------------------------------------------------------------------
 * subroutine:  fft
 * multivariate complex fourier transform, computed in place
 * using mixed-radix fast fourier transform algorithm.
 *-----------------------------------------------------------------------
 *
 *      this is the call from C:
 *              fft_(anal,banal,&one,&N2,&one,&mtwo);
 *      CHANGED TO:-
 *              fft_(anal,banal,one,N2,one,mtwo);
 */

void
fft_(double *a, double *b, int nseg, int n, int nspn, int isn)
  /*    *a,       pointer to array 'anal'  */
  /*    *b;       pointer to array 'banal' */
{
    int nfac[/*16*/32];               /*  (2023)WAS 16. These are one bigger than needed   */
                                /*  because wish to use Fortran array  */
                                /* index which runs 1 to n, not 0 to n */

    int         m = 0,
                nf,
                k,
                kt,
                ntot,
                j,
                jj,
                maxf, maxp=-1;

/* work space pointers */
    double       *at, *ck, *bt, *sk;
    int *np;


/* reduce the pointers to input arrays - by doing this, FFT uses FORTRAN
   indexing but retains compatibility with C arrays */
    a--;        b--;

/*
 * determine the factors of n
 */
    k=nf=abs(n);
    if (nf==1)
      return;

    nspn=abs(nf*nspn);
    ntot=abs(nspn*nseg);

    if (isn*ntot == 0) {
      printf("\nerror - zero in fft parameters %d %d %d %d",
              nseg, n, nspn, isn);
      return;
    }
    for (m=0; !(k%16); nfac[++m]=4,k/=16);
    for (j=3,jj=9; jj<=k; j+=2,jj=j*j)
      for (; !(k%jj); nfac[++m]=j,k/=jj);

    if (k<=4) {
      kt = m;
      nfac[m+1] = k;
      if (k != 1)
        m++;
    }
    else {
      if (k%4==0) {
        nfac[++m]=2;
        k/=4;
      }

      kt = m;
      maxp = (kt+kt+2 > k-1 ? kt+kt+2 : k-1);
      for (j=2; j<=k; j=1+((j+1)/2)*2)
        if (k%j==0) {
          nfac[++m]=j;
          k/=j;
        }
    }
    if (m <= kt+1)
      maxp = m + kt + 1;
    if (m+kt > 15) {
      printf("\nerror - fft parameter n has more than 15 factors : %d", n);
      return;
    }
    if (kt!=0) {
      j = kt;
      while (j)
        nfac[++m]=nfac[j--];
    }
    maxf = nfac[m-kt];
    if (kt > 0 && maxf <nfac[kt])
      maxf = nfac[kt];

/*  allocate workspace - assume no errors! */
    at = (double *) calloc(maxf,sizeof(double));
    ck = (double *) calloc(maxf,sizeof(double));
    bt = (double *) calloc(maxf,sizeof(double));
    sk = (double *) calloc(maxf,sizeof(double));
    np = (int *) calloc(maxp,sizeof(int));

/* decrement pointers to allow FORTRAN type usage in fftmx */
    at--;       bt--;   ck--;   sk--;   np--;

/* call fft driver */

    fftmx(a,b,ntot,nf,nspn,isn,m,&kt,at,ck,bt,sk,np,nfac);

/* restore pointers before releasing */
    at++;       bt++;   ck++;   sk++;   np++;

/* release working storage before returning - assume no problems */
    free(at);
    free(sk);
    free(bt);
    free(ck);
    free(np);
    return;
}

/*
 *-----------------------------------------------------------------------
 * subroutine:  fftmx
 * called by subroutine 'fft' to compute mixed-radix fourier transform
 *-----------------------------------------------------------------------
 */
void
fftmx(double *a, double *b, int ntot, int n, int nspan, int isn, int m,
      int *kt, double *at, double *ck, double *bt, double *sk, int *np, int nfac[])
{
    int i,inc,
      j,jc,jf, jj,
      k, k1, k2, k3 = 0, k4,
      kk,klim,ks,kspan, kspnn,
      lim,
      maxf,mm,
      nn,nt;
    double  aa, aj, ajm, ajp, ak, akm, akp,
      bb, bj, bjm, bjp, bk, bkm, bkp,
      c1, c2 = 0.0, c3 = 0.0, c72, cd,
      dr,
      rad,
      sd, s1, s2=0.0, s3=0.0, s72, s120;

    double      xx;     /****** ADDED APRIL 1991 *********/
    inc=abs(isn);
    nt = inc*ntot;
    ks = inc*nspan;
/******************* REPLACED MARCH 29: ***********************
                                        rad = atan((double)1.0);
**************************************************************/
    rad = 0.785398163397448278900;
/******************* REPLACED MARCH 29: ***********************
                                        s72 = rad/0.625;
                                        c72 = cos(s72);
                                        s72 = sin(s72);
**************************************************************/
    c72 = 0.309016994374947451270;
    s72 = 0.951056516295153531190;
/******************* REPLACED MARCH 29: ***********************
                                        s120 = sqrt((double)0.75);
**************************************************************/
    s120 = 0.866025403784438707600;

/* scale by 1/n for isn > 0 ( reverse transform ) */

    if (isn < 0) {
      s72 = -s72;
      s120 = -s120;
      rad = -rad;}
    else {
      ak = 1.0/(double)n;
      for (j=1; j<=nt;j += inc) {
        a[j] *= (double)ak;
        b[j] *= (double)ak;
      }
    }
    kspan = ks;
    nn = nt - inc;
    jc = ks/n;

/* sin, cos values are re-initialised each lim steps  */

    lim = 32;
    klim = lim * jc;
    i = 0;
    jf = 0;
    maxf = m - (*kt);
    maxf = nfac[maxf];
    if ((*kt) > 0 && maxf < nfac[*kt])
      maxf = nfac[*kt];

/*
 * compute fourier transform
 */

 lbl40:
    dr = (8.0 * (double)jc)/((double)kspan);
/*************************** APRIL 1991 POW & POW2 not WORKING.. REPLACE *******
                    cd = 2.0 * (pow2 ( sin((double)0.5 * dr * rad)) );
*******************************************************************************/
    xx =  sin((double)0.5 * dr * rad);
    cd = 2.0 * xx * xx;
    sd = sin(dr * rad);
    kk = 1;
    if (nfac[++i]!=2) goto lbl110;
/*
 * transform for factor of 2 (including rotation factor)
 */
    kspan /= 2;
    k1 = kspan + 2;
    do {
      do {
        k2 = kk + kspan;
        ak = a[k2];
        bk = b[k2];
        a[k2] = (a[kk]) - (double)ak;
        b[k2] = (b[kk]) - (double)bk;
        a[kk] = (a[kk]) + (double)ak;
        b[kk] = (b[kk]) + (double)bk;
        kk = k2 + kspan;
      } while (kk <= nn);
      kk -= nn;
    } while (kk <= jc);
    if (kk > kspan) goto lbl350;
 lbl60:
    c1 = 1.0 - cd;
    s1 = sd;
    mm = (k1/2 < klim ? k1/2 :klim);
    goto lbl80;
 lbl70:
    ak = c1 - ((cd*c1)+(sd*s1));
    s1 = ((sd*c1)-(cd*s1)) + s1;
    c1 = ak;
lbl80:
    do {
      do {
        k2 = kk + kspan;
        ak = a[kk] - a[k2];
        bk = b[kk] - b[k2];
        a[kk] = a[kk] + a[k2];
        b[kk] = b[kk] + b[k2];
        a[k2] = (double)((c1 * ak) - (s1 * bk));
        b[k2] = (double)((s1 * ak) + (c1 * bk));
        kk = k2 + kspan;
      } while (kk < nt);
      k2 = kk - nt;
      c1 = -c1;
      kk = k1 - k2;
    } while (kk > k2);
    kk += jc;
    if (kk <= mm) goto lbl70;
    if (kk < k2)  goto lbl90;
    k1 += (inc + inc);
    kk = ((k1-kspan)/2) + jc;
    if (kk <= (jc+jc)) goto lbl60;
    goto lbl40;
 lbl90:
    s1 = ((double)((kk-1)/jc)) * dr * rad;
    c1 = cos(s1);
    s1 = sin(s1);
    mm = (k1/2 < mm+klim ? k1/2 : mm+klim);
    goto lbl80;
/*
 * transform for factor of 3 (optional code)
 */


 lbl100:
    k1 = kk + kspan;
    k2 = k1 + kspan;
    ak = a[kk];
    bk = b[kk];
    aj = a[k1] + a[k2];
    bj = b[k1] + b[k2];
    a[kk] = (double)(ak + aj);
    b[kk] = (double)(bk + bj);
    ak += (-0.5 * aj);
    bk += (-0.5 * bj);
    aj = (a[k1] - a[k2]) * s120;
    bj = (b[k1] - b[k2]) * s120;
    a[k1] = (double)(ak - bj);
    b[k1] = (double)(bk + aj);
    a[k2] = (double)(ak + bj);
    b[k2] = (double)(bk - aj);
    kk = k2 + kspan;
    if (kk < nn)     goto lbl100;
    kk -= nn;
    if (kk <= kspan) goto lbl100;
    goto lbl290;

/*
 * transform for factor of 4
 */

 lbl110:
    if (nfac[i] != 4) goto lbl230;
    kspnn = kspan;
    kspan = kspan/4;
 lbl120:
    c1 = 1.0;
    s1 = 0;
    mm = (kspan < klim ? kspan : klim);
    goto lbl150;
 lbl130:
    c2 = c1 - ((cd*c1)+(sd*s1));
    s1 = ((sd*c1)-(cd*s1)) + s1;
/*
 * the following three statements compensate for truncation
 * error.  if rounded arithmetic is used, substitute
 * c1=c2
 *
 * c1 = (0.5/(pow2(c2)+pow2(s1))) + 0.5;
 * s1 = c1*s1;
 * c1 = c1*c2;
 */
    c1 = c2;
 lbl140:
    c2 = (c1 * c1) - (s1 * s1);
    s2 = c1 * s1 * 2.0;
    c3 = (c2 * c1) - (s2 * s1);
    s3 = (c2 * s1) + (s2 * c1);
 lbl150:
    k1 = kk + kspan;
    k2 = k1 + kspan;
    k3 = k2 + kspan;
    akp = a[kk] + a[k2];
    akm = a[kk] - a[k2];
    ajp = a[k1] + a[k3];
    ajm = a[k1] - a[k3];
    a[kk] = (double)(akp + ajp);
    ajp = akp - ajp;
    bkp = b[kk] + b[k2];
    bkm = b[kk] - b[k2];
    bjp = b[k1] + b[k3];
    bjm = b[k1] - b[k3];
    b[kk] = (double)(bkp + bjp);
    bjp = bkp - bjp;
    if (isn < 0) goto lbl180;
    akp = akm - bjm;
    akm = akm + bjm;
    bkp = bkm + ajm;
    bkm = bkm - ajm;
    if (s1 == 0.0) goto lbl190;
 lbl160:
    a[k1] = (double)((akp*c1) - (bkp*s1));
    b[k1] = (double)((akp*s1) + (bkp*c1));
    a[k2] = (double)((ajp*c2) - (bjp*s2));
    b[k2] = (double)((ajp*s2) + (bjp*c2));
    a[k3] = (double)((akm*c3) - (bkm*s3));
    b[k3] = (double)((akm*s3) + (bkm*c3));
    kk = k3 + kspan;
    if (kk <= nt)   goto lbl150;
 lbl170:
    kk -= (nt - jc);
    if (kk <= mm)   goto lbl130;
    if (kk < kspan) goto lbl200;
    kk -= (kspan - inc);
    if (kk <= jc)   goto lbl120;
    if (kspan==jc)  goto lbl350;
    goto lbl40;
lbl180:
    akp = akm + bjm;
    akm = akm - bjm;
    bkp = bkm - ajm;
    bkm = bkm + ajm;
    if (s1 != 0.0)  goto lbl160;
 lbl190:
    a[k1] = (double)akp;
    b[k1] = (double)bkp;
    a[k2] = (double)ajp;
    b[k2] = (double)bjp;
    a[k3] = (double)akm;
    b[k3] = (double)bkm;
    kk = k3 + kspan;
    if (kk <= nt) goto lbl150;
    goto lbl170;
 lbl200:
    s1 = ((double)((kk-1)/jc)) * dr * rad;
    c1 = cos(s1);
    s1 = sin(s1);
    mm = (kspan < mm+klim ? kspan : mm+klim);
    goto lbl140;

/*
 * transform for factor of 5 (optional code)
 */

 lbl210:
    c2 = (c72*c72) - (s72*s72);
    s2 = 2.0 * c72 * s72;
 lbl220:
    k1 = kk + kspan;
    k2 = k1 + kspan;
    k3 = k2 + kspan;
    k4 = k3 + kspan;
    akp = a[k1] + a[k4];
    akm = a[k1] - a[k4];
    bkp = b[k1] + b[k4];
    bkm = b[k1] - b[k4];
    ajp = a[k2] + a[k3];
    ajm = a[k2] - a[k3];
    bjp = b[k2] + b[k3];
    bjm = b[k2] - b[k3];
    aa = a[kk];
    bb = b[kk];
    a[kk] = (double)(aa + akp + ajp);
    b[kk] = (double)(bb + bkp + bjp);
    ak = (akp*c72) + (ajp*c2) + aa;
    bk = (bkp*c72) + (bjp*c2) + bb;
    aj = (akm*s72) + (ajm*s2);
    bj = (bkm*s72) + (bjm*s2);
    a[k1] = (double)(ak - bj);
    a[k4] = (double)(ak + bj);
    b[k1] = (double)(bk + aj);
    b[k4] = (double)(bk - aj);
    ak = (akp*c2) + (ajp*c72) + aa;
    bk = (bkp*c2) + (bjp*c72) + bb;
    aj = (akm*s2) - (ajm*s72);
    bj = (bkm*s2) - (bjm*s72);
    a[k2] = (double)(ak - bj);
    a[k3] = (double)(ak + bj);
    b[k2] = (double)(bk + aj);
    b[k3] = (double)(bk - aj);
    kk = k4 + kspan;
    if (kk < nn)     goto lbl220;
    kk -= nn;
    if (kk <= kspan) goto lbl220;
    goto lbl290;

/*
 * transform for odd factors
 */

 lbl230:
    k = nfac[i];
    kspnn = kspan;
    kspan /= k;
    if (k==3)   goto lbl100;
    if (k==5)   goto lbl210;
    if (k==jf)  goto lbl250;
    jf = k;
    s1 = rad/(((double)(k))/8.0);
    c1 = cos(s1);
    s1 = sin(s1);
    ck[jf] = 1.0;
    sk[jf] = 0.0;
    for (j=1; j<k ; j++) {
      ck[j] = (double)((ck[k])*c1 + (sk[k])*s1);
      sk[j] = (double)((ck[k])*s1 - (sk[k])*c1);
      k--;
      ck[k] = ck[j];
      sk[k] = -(sk[j]);
    }
 lbl250:
    k1 = kk;
    k2 = kk + kspnn;
    aa = a[kk];
    bb = b[kk];
    ak = aa;
    bk = bb;
    j = 1;
    k1 += kspan;
    do {
      k2 -= kspan;
      j++;
      at[j] = a[k1] + a[k2];
      ak = at[j] + ak;
      bt[j] = b[k1] + b[k2];
      bk = bt[j] + bk;
      j++;
      at[j] = a[k1] - a[k2];
      bt[j] = b[k1] - b[k2];
      k1 += kspan;
    } while (k1 < k2);
    a[kk] = (double)ak;
    b[kk] = (double)bk;
    k1 = kk;
    k2 = kk + kspnn;
    j = 1;
 lbl270:
    k1 += kspan;
    k2 -= kspan;
    jj = j;
    ak = aa;
    bk = bb;
    aj = 0.0;
    bj = 0.0;
    k = 1;
    do {
      k++;
      ak = (at[k] * ck[jj]) + ak;
      bk = (bt[k] * ck[jj]) + bk;
      k++;
      aj = (at[k] * sk[jj]) + aj;
      bj = (bt[k] * sk[jj]) + bj;
      jj += j;
      if (jj > jf)
        jj -= jf;
    } while (k < jf);
    k = jf - j;
    a[k1] = (double)(ak - bj);
    b[k1] = (double)(bk + aj);
    a[k2] = (double)(ak + bj);
    b[k2] = (double)(bk - aj);
    j++;
    if (j < k)     goto lbl270;
    kk += kspnn;
    if (kk <= nn)  goto lbl250;
    kk -= nn;
    if (kk<=kspan) goto lbl250;

/*
 * multiply by rotation factor (except for factors of 2 and 4)
 */

 lbl290:
    if (i==m) goto lbl350;
    kk = jc + 1;
 lbl300:
    c2 = 1.0 - cd;
    s1 = sd;
    mm = (kspan < klim ? kspan : klim);
    goto lbl320;
 lbl310:
    c2 = c1 - ((cd*c1) + (sd*s1));
    s1 = s1 + ((sd*c1) - (cd*s1));
 lbl320:
    c1 = c2;
    s2 = s1;
    kk += kspan;
 lbl330:
    ak = a[kk];
    a[kk] = (double)((c2*ak) - (s2 * b[kk]));
    b[kk] = (double)((s2*ak) + (c2 * b[kk]));
    kk += kspnn;
    if (kk <= nt) goto lbl330;
    ak = s1*s2;
    s2 = (s1*c2) + (c1*s2);
    c2 = (c1*c2) - ak;
    kk -= (nt - kspan);
    if (kk <= kspnn) goto lbl330;
    kk -= (kspnn - jc);
    if (kk <= mm)   goto lbl310;
    if (kk < kspan) goto lbl340;
    kk -= (kspan - jc - inc);
    if (kk <= (jc+jc)) goto lbl300;
    goto lbl40;
 lbl340:
    s1 = ((double)((kk-1)/jc)) * dr * rad;
    c2 = cos(s1);
    s1 = sin(s1);
    mm = (kspan < mm+klim ?  kspan :mm+klim);
    goto lbl320;

/*
 * permute the results to normal order---done in two stages
 * permutation for square factors of n
 */

 lbl350:
    np[1] = ks;
    if (!(*kt)) goto lbl440;
    k = *kt + *kt + 1;
    if (m < k)
      k--;
    np[k+1] = jc;
    for (j=1; j < k; j++,k--) {
      np[j+1] = np[j] / nfac[j];
      np[k] = np[k+1] * nfac[j];
    }
    k3 = np[k+1];
    kspan = np[2];
    kk = jc + 1;
    k2 = kspan + 1;
    j = 1;
    if (n != ntot) goto lbl400;
/*
 * permutation for single-variate transform (optional code)
 */
 lbl370:
    do {
      ak = a[kk];
      a[kk] = a[k2];
      a[k2] = (double)ak;
      bk = b[kk];
      b[kk] = b[k2];
      b[k2] = (double)bk;
      kk += inc;
      k2 += kspan;
    } while (k2 < ks);
lbl380:
    do {
      k2 -= np[j++];
      k2 += np[j+1];
    } while (k2 > np[j]);
    j = 1;
 lbl390:
    if (kk < k2) {
      goto lbl370;
    }
    kk += inc;
    k2 += kspan;
    if (k2 < ks) goto lbl390;
    if (kk < ks) goto lbl380;
    jc = k3;
    goto lbl440;
/*
 * permutation for multivariate transform
 */
 lbl400:
    do {
      do {
        k = kk + jc;
        do {
          ak = a[kk];
          a[kk] = a[k2];
          a[k2] = (double)ak;
          bk = b[kk];
          b[kk] = b[k2];
          b[k2] = (double)bk;
          kk += inc;
          k2 += inc;
        } while (kk < k);
        kk += (ks - jc);
        k2 += (ks - jc);
      } while (kk < nt);
      k2 -= (nt - kspan);
      kk -= (nt - jc);
    } while (k2 < ks);
 lbl420:
    do {
      k2 -= np[j++];
      k2 += np[j+1];
    } while (k2 > np[j]);
    j = 1;
 lbl430:
    if (kk < k2)         goto lbl400;
    kk += jc;
    k2 += kspan;
    if (k2 < ks)      goto lbl430;
    if (kk < ks)      goto lbl420;
    jc = k3;
 lbl440:
    if ((2*(*kt))+1 >= m)
      return;

    kspnn = *(np + *(kt) + 1);
    j = m - *kt;
    nfac[j+1] = 1;
 lbl450:
    nfac[j] = nfac[j] * nfac[j+1];
    j--;
    if (j != *kt) goto lbl450;
    *kt = *(kt) + 1;
    nn = nfac[*kt] - 1;
    jj = 0;
    j = 0;
    goto lbl480;
 lbl460:
    jj -= k2;
    k2 = kk;
    kk = nfac[++k];
 lbl470:
    jj += kk;
    if (jj >= k2) goto lbl460;
    np[j] = jj;
 lbl480:
    k2 = nfac[*kt];
    k = *kt + 1;
    kk = nfac[k];
    j++;
    if (j <= nn) goto lbl470;
/* Determine permutation cycles of length greater than 1 */
    j = 0;
    goto lbl500;
 lbl490:
    k = kk;
    kk = np[k];
    np[k] = -kk;
    if (kk != j) goto lbl490;
    k3 = kk;
 lbl500:
    kk = np[++j];
    if (kk < 0)  goto lbl500;
    if (kk != j) goto lbl490;
    np[j] = -j;
    if (j != nn) goto lbl500;
    maxf *= inc;
    /* Perform reordering following permutation cycles */
    goto lbl570;
 lbl510:
    j--;
    if (np[j] < 0) goto lbl510;
    jj = jc;
 lbl520:
    kspan = jj;
    if (jj > maxf)
      kspan = maxf;
    jj -= kspan;
    k = np[j];
    kk = (jc*k) + i + jj;
    k1 = kk + kspan;
    k2 = 0;
 lbl530:
    k2++;
    at[k2] = a[k1];
    bt[k2] = b[k1];
    k1 -= inc;
    if (k1 != kk) goto lbl530;
 lbl540:
    k1 = kk + kspan;
    k2 = k1 - (jc * (k + np[k]));
    k = -(np[k]);
 lbl550:
    a[k1] = a[k2];
    b[k1] = b[k2];
    k1 -= inc;
    k2 -= inc;
    if (k1 != kk) goto lbl550;
    kk = k2;
    if (k != j)   goto lbl540;
    k1 = kk + kspan;
    k2 = 0;
 lbl560:
    k2++;
    a[k1] = at[k2];
    b[k1] = bt[k2];
    k1 -= inc;
    if (k1 != kk) goto lbl560;
    if (jj)       goto lbl520;
    if (j  != 1)  goto lbl510;
lbl570:
    j = k3 + 1;
    nt -= kspnn;
    i = nt - inc + 1;
    if (nt >= 0)  goto lbl510;
    return;
}


/*
 *-----------------------------------------------------------------------
 * subroutine:
  reals
 * used with 'fft' to compute fourier transform or inverse for real data
 *-----------------------------------------------------------------------
 *      this is the call from C:

 *              reals_(anal,banal,N2,mtwo);
 *      which has been changed from CARL call
 *              reals_(anal,banal,&N2,&mtwo);
 */

void
reals_(double *a, double *b, int n, int isn)

  /*    *a,       a refers to an array of floats 'anal'   */
  /*    *b;       b refers to an array of floats 'banal'  */
/* See IEEE book for a long comment here on usage */

{
    int inc,
      j,
      k,
      lim,
      mm,ml,
      nf,nk,nh;

    double      aa,ab,
      ba,bb,
      cd,cn,
      dr,
      em,
      rad,re,
      sd,sn;
    double      xx;     /******* ADDED APRIL 1991 ******/
    /* adjust  input array pointers (called from C) */
    a--;        b--;
    inc=abs(isn);
    nf=abs(n);
    if (nf*isn==0) {
      printf("\nerror - zero in reals parameters : %d : %d ",n,isn);
      return;
    }
    nk = (nf*inc) + 2;
    nh = nk/2;
/*****************************
        rad  = atan((double)1.0);
******************************/
    rad = 0.785398163397448278900;
    dr = -4.0/(double)(nf);
/********************************** POW2 REMOVED APRIL 1991 *****************
                                cd = 2.0 * (pow2(sin((double)0.5 * dr * rad)));
*****************************************************************************/
    xx = sin((double)0.5 * dr * rad);
    cd = 2.0 * xx * xx;
    sd = sin(dr * rad);
/*
 * sin,cos values are re-initialised each lim steps
 */
    lim = 32;
    mm = lim;
    ml = 0;
    sn = 0.0;
    if (isn<0) {
      cn = 1.0;
      a[nk-1] = a[1];
      b[nk-1] = b[1]; }
    else {
      cn = -1.0;
      sd = -sd;
    }
    for (j=1;j<=nh;j+=inc)      {
      k = nk - j;
      aa = a[j] + a[k];
      ab = a[j] - a[k];
      ba = b[j] + b[k];
      bb = b[j] - b[k];
      re = (cn*ba) + (sn*ab);
      em = (sn*ba) - (cn*ab);
      b[k] = (double)((em-bb)*0.5);
      b[j] = (double)((em+bb)*0.5);
      a[k] = (double)((aa-re)*0.5);
      a[j] = (double)((aa+re)*0.5);
      ml++;
      if (ml!=mm) {
        aa = cn - ((cd*cn)+(sd*sn));
        sn = ((sd*cn) - (cd*sn)) + sn;
        cn = aa;}
      else {
        mm +=lim;
        sn = ((double)ml) * dr * rad;
        cn = cos(sn);
        if (isn>0)
          cn = -cn;
        sn = sin(sn);
      }
    }
    return;
}