dpfafft.c

该程序是用vc开发的对动态数组进行管理的DLL

				t26r = t21r-t22r;
				t26i = t21i-t22i;
				t27r = t9r+t24r;
				t27i = t9i+t24i;
				t28r = t10r+t23r;
				t28i = t10i+t23i;
				t29r = t9r-t24r;
				t29i = t9i-t24i;
				t30r = t10r-t23r;
				t30i = t10i-t23i;
				t31r = t5r+t17r;
				t31i = t5i+t17i;
				t32r = t11r+t25r;
				t32i = t11i+t25i;
				t33r = t3r+t18r;
				t33i = t3i+t18i;
				t34r = c2*t29r-c6*t30r;
				t34i = c2*t29i-c6*t30i;
				t35r = t3r-t18r;
				t35i = t3i-t18i;
				t36r = c7*t27r-c3*t28r;
				t36i = c7*t27i-c3*t28i;
				t37r = t4r+t19r;
				t37i = t4i+t19i;
				t38r = c3*t27r+c7*t28r;
				t38i = c3*t27i+c7*t28i;
				t39r = t4r-t19r;
				t39i = t4i-t19i;
				t40r = c6*t29r+c2*t30r;
				t40i = c6*t29i+c2*t30i;
				t41r = c4*(t12r-t26r);
				t41i = c4*(t12i-t26i);
				t42r = c5*(t12r+t26r);
				t42i = c5*(t12i+t26i);
				y1r = t33r+t34r;
				y1i = t33i+t34i;
				y2r = t6r+t41r;
				y2i = t6i+t41i;
				y3r = t35r+t40r;
				y3i = t35i+t40i;
				y4r = t5r-t17r;
				y4i = t5i-t17i;
				y5r = t35r-t40r;
				y5i = t35i-t40i;
				y6r = t6r-t41r;
				y6i = t6i-t41i;
				y7r = t33r-t34r;
				y7i = t33i-t34i;
				y9r = t38r-t37r;
				y9i = t38i-t37i;
				y10r = t42r-t20r;
				y10i = t42i-t20i;
				y11r = t36r+t39r;
				y11i = t36i+t39i;
				y12r = c1*(t11r-t25r);
				y12i = c1*(t11i-t25i);
				y13r = t36r-t39r;
				y13i = t36i-t39i;
				y14r = t42r+t20r;
				y14i = t42i+t20i;
				y15r = t38r+t37r;
				y15i = t38i+t37i;
				z[j00] = t31r+t32r;
				z[j00+1] = t31i+t32i;
				z[j01] = y1r-y15i;
				z[j01+1] = y1i+y15r;
				z[j2] = y2r-y14i;
				z[j2+1] = y2i+y14r;
				z[j3] = y3r-y13i;
				z[j3+1] = y3i+y13r;
				z[j4] = y4r-y12i;
				z[j4+1] = y4i+y12r;
				z[j5] = y5r-y11i;
				z[j5+1] = y5i+y11r;
				z[j6] = y6r-y10i;
				z[j6+1] = y6i+y10r;
				z[j7] = y7r-y9i;
				z[j7+1] = y7i+y9r;
				z[j8] = t31r-t32r;
				z[j8+1] = t31i-t32i;
				z[j9] = y7r+y9i;
				z[j9+1] = y7i-y9r;
				z[j10] = y6r+y10i;
				z[j10+1] = y6i-y10r;
				z[j11] = y5r+y11i;
				z[j11+1] = y5i-y11r;
				z[j12] = y4r+y12i;
				z[j12+1] = y4i-y12r;
				z[j13] = y3r+y13i;
				z[j13+1] = y3i-y13r;
				z[j14] = y2r+y14i;
				z[j14+1] = y2i-y14r;
				z[j15] = y1r+y15i;
				z[j15+1] = y1i-y15r;
				jt = j15+2;
				j15 = j14+2;
				j14 = j13+2;
				j13 = j12+2;
				j12 = j11+2;
				j11 = j10+2;
				j10 = j9+2;
				j9 = j8+2;
				j8 = j7+2;
				j7 = j6+2;
				j6 = j5+2;
				j5 = j4+2;
				j4 = j3+2;
				j3 = j2+2;
				j2 = j01+2;
				j01 = j00+2;
				j00 = jt;
			}
			continue;
		}
	}
}

void pfacr_d (int isign, int n, real_complex cz[], real rz[])
/*****************************************************************************
Prime factor fft:  real_complex to real transform
******************************************************************************
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform (see notes below)
cz          array[n/2+1] of real_complex values (may be equivalenced to rz)

Output:
rz          array[n] of real values (may be equivalenced to cz)
******************************************************************************
Notes:
Because pfacr uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.
******************************************************************************
References:  
Press et al, 1988, Numerical Recipes in C, p. 417.

Also, see notes and references for function pfacc.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/13/89
*****************************************************************************/
{
    int i,ir,ii,jr,ji,no2;
    real *z,tempr,tempi,sumr,sumi,difr,difi;
    double wr,wi,wpr,wpi,wtemp,theta;

    /* copy input to output and fix dc and nyquist */
    z = (real*)cz;
    for (i=2; i<n; i++)
        rz[i] = z[i];
    rz[1] = z[0]-z[n];
    rz[0] = z[0]+z[n];
    z = rz;

    /* initialize cosine-sine recurrence */
    theta = 2.0*PI/(double)n;
    if (isign>0) theta = -theta;
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0+wpr;
    wi = wpi;

    /* twiddle */
    no2 = n/2;
    for (ir=2,ii=3,jr=n-2,ji=n-1; ir<=no2; ir+=2,ii+=2,jr-=2,ji-=2) {
        sumr = z[ir]+z[jr];
        sumi = z[ii]+z[ji];
        difr = z[ir]-z[jr];
        difi = z[ii]-z[ji];
        tempr = wi*difr-wr*sumi;
        tempi = wi*sumi+wr*difr;
        z[ir] = sumr+tempr;
        z[ii] = difi+tempi;
        z[jr] = sumr-tempr;
        z[ji] = tempi-difi;
        wtemp = wr;
        wr += wr*wpr-wi*wpi;
        wi += wi*wpr+wtemp*wpi;
    }

    /* do real_complex to real_complex transform */
    pfacc_d(isign,n/2,(real_complex*)z);
}

void pfarc_d (int isign, int n, real rz[], real_complex cz[])
/*****************************************************************************
Prime factor fft:  real to real_complex transform
******************************************************************************
Input:
isign       sign of isign is the sign of exponent in fourier kernel
n           length of transform; must be even (see notes below)
rz          array[n] of real values (may be equivalenced to cz)

Output:
cz          array[n/2+1] of real_complex values (may be equivalenced to rz)
******************************************************************************
Notes:
Because pfarc uses pfacc to do most of the work, n must be even 
and n/2 must be a valid length for pfacc.  The simplest way to
obtain a valid n is via n = npfar(nmin).  A more optimal n can be 
obtained with npfaro.
******************************************************************************
References:  
Press et al, 1988, Numerical Recipes in C, p. 417.

Also, see notes and references for function pfacc.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/13/89
*****************************************************************************/
{
    int i,ir,ii,jr,ji,no2;
    real *z,tempr,tempi,sumr,sumi,difr,difi;
    double wr,wi,wpr,wpi,wtemp,theta;

    /* copy input to output while scaling */
    z = (real*)cz;
    for (i=0; i<n; i++)
        z[i] = 0.5*rz[i];

    /* do real_complex to real_complex transform */
    pfacc_d(isign,n/2,cz);

    /* fix dc and nyquist */
    z[n] = 2.0*(z[0]-z[1]);
    z[0] = 2.0*(z[0]+z[1]);
    z[n+1] = 0.0;
    z[1] = 0.0;

    /* initialize cosine-sine recurrence */
    theta = 2.0*PI/(double)n;
    if (isign<0) theta = -theta;
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0+wpr;
    wi = wpi;

    /* twiddle */
    no2 = n/2;
    for (ir=2,ii=3,jr=n-2,ji=n-1; ir<=no2; ir+=2,ii+=2,jr-=2,ji-=2) {
        sumr = z[ir]+z[jr];
        sumi = z[ii]+z[ji];
        difr = z[ir]-z[jr];
        difi = z[ii]-z[ji];
        tempr = wi*difr+wr*sumi;
        tempi = wi*sumi-wr*difr;
        z[ir] = sumr+tempr;
        z[ii] = difi+tempi;
        z[jr] = sumr-tempr;
        z[ji] = tempi-difi;
        wtemp = wr;
        wr += wr*wpr-wi*wpi;
        wi += wi*wpr+wtemp*wpi;
    }
}

void pfamcc_d (int isign, int n, int nt, int k, int kt, real_complex cz[])
/*****************************************************************************
Prime factor fft:  multiple real_complex to real_complex transforms, in place
******************************************************************************
Input:
isign       	sign of isign is the sign of exponent in fourier kernel
n           	number of real_complex elements per transform (see notes below)
nt          	number of transforms
k           	stride in real_complex elements within transforms
kt          	stride in real_complex elements between transforms
z           	array of real_complex elements to be transformed in place

Output:
z		array of real_complex elements transformed
******************************************************************************
Notes:
n must be factorable into mutually prime factors taken 
from the set {2,3,4,5,7,8,9,11,13,16}.  in other words,
    n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
where
    0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
is required for pfamcc to yield meaningful results.  this
restriction implies that n is restricted to the range
    1 <= n <= 720720 (= 5*7*9*11*13*16)

To perform a two-dimensional transform of an n1 by n2 real_complex array 
(assuming that both n1 and n2 are valid "n"), stored with n1 fast 
and n2 slow:
    pfamcc(isign,n1,n2,1,n1,z); (to transform 1st dimension)
    pfamcc(isign,n2,n1,n1,1,z); (to transform 2nd dimension)
******************************************************************************
References:  
Temperton, C., 1985, Implementation of a self-sorting
in-place prime factor fft algorithm:  Journal of
Computational Physics, v. 58, p. 283-299.

Temperton, C., 1988, A new set of minimum-add rotated
rotated dft modules: Journal of Computational Physics,
v. 75, p. 190-198.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/15/89
*****************************************************************************/
{
    static int kfax[] = { 16,13,11,9,8,7,5,4,3,2 };
    register real *z=(real*)cz;
    register int j00,j01,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14,j15;
    int nleft,jfax,ifac,jfac,iinc,imax,ndiv,m,mm=0,mu=0,l,istep,jstep,
        jt,i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15,it;
    real t1r,t1i,t2r,t2i,t3r,t3i,t4r,t4i,t5r,t5i,
        t6r,t6i,t7r,t7i,t8r,t8i,t9r,t9i,t10r,t10i,
        t11r,t11i,t12r,t12i,t13r,t13i,t14r,t14i,t15r,t15i,
        t16r,t16i,t17r,t17i,t18r,t18i,t19r,t19i,t20r,t20i,
        t21r,t21i,t22r,t22i,t23r,t23i,t24r,t24i,t25r,t25i,
        t26r,t26i,t27r,t27i,t28r,t28i,t29r,t29i,t30r,t30i,
        t31r,t31i,t32r,t32i,t33r,t33i,t34r,t34i,t35r,t35i,
        t36r,t36i,t37r,t37i,t38r,t38i,t39r,t39i,t40r,t40i,
        t41r,t41i,t42r,t42i,
        y1r,y1i,y2r,y2i,y3r,y3i,y4r,y4i,y5r,y5i,
        y6r,y6i,y7r,y7i,y8r,y8i,y9r,y9i,y10r,y10i,
        y11r,y11i,y12r,y12i,y13r,y13i,y14r,y14i,y15r,y15i,
        c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12;

    /* determine step within and between transforms */
    istep = 2*k;
    jstep = 2*kt;

    /* keep track of n left after dividing by factors */
    nleft = n;

    /* begin loop over possible factors (from biggest to smallest) */
    for (jfax=0; jfax<NFAX; jfax++) {

        /* skip if not a mutually prime factor of n */
        ifac = kfax[jfax];
        ndiv = nleft/ifac;
        if (ndiv*ifac!=nleft) continue;
 
        /* update n left and determine n divided by factor */
        nleft = ndiv;
        m = n/ifac;
 
        /* determine rotation factor mu and stride mm */
        for (jfac=1; jfac<=ifac; jfac++) {
            mu = jfac;
            mm = jfac*m;
            if (mm%ifac==1) break;
        }
 
        /* adjust rotation factor for sign of transform */
        if (isign<0) mu = ifac-mu;
 
        /* compute stride, limit, and pointers */
        iinc = istep*mm;
        imax = istep*n;
        i0 = 0;
        i1 = i0+iinc;

        /* if factor is 2 */
        if (ifac==2) {
            for (l=0; l<m; l++) {
                j00 = i0;
                j01 = i1;
                for (jt=0; jt<nt; jt++) {
                    t1r = z[j00]-z[j01];
                    t1i = z[j00+1]-z[j01+1];
                    z[j00] = z[j00]+z[j01];
                    z[j00+1] = z[j00+1]+z[j01+1];
                    z[j01] = t1r;
                    z[j01+1] = t1i;
                    j00 += jstep;
                    j01 += jstep;
                }
                it = i1+istep;
                i1 = i0+istep;
                i0 = it;
            }
            continue;
        }
        i2 = i1+iinc;
        if (i2>=imax) i2 = i2-imax;

        /* if factor is 3 */
        if (ifac==3) {
            if (mu==1)
                c1 = P866;
            else
                c1 = -P866;
            for (l=0; l<m; l++) {
                j00 = i0;
                j01 = i1;
                j2 = i2;
                for (jt=0; jt<nt; jt++) {
                    t1r = z[j01]+z[j2];
                    t1i = z[j01+1]+z[j2+1];
                    y1r = z[j00]-0.5*t1r;
                    y1i = z[j00+1]-0.5*t1i;
                    y2r = c1*(z[j01]-z[j2]);
                    y2i = c1*(z[j01+1]-z[j2+1]);
                    z[j00] = z[j00]+t1r;
                    z[j00+1] = z[j00+1]+t1i;
                    z[j01] = y1r-y2i;
                    z[j01+1] = y1i+y2r;
                    z[j2] = y1r+y2i;
                    z[j2+1] = y1i-y2r;
                    j00 += jstep;
                    j01 += jstep;
                    j2 += jstep;
                }
                it = i2+istep;
                i2 = i1+istep;
                i1 = i0+istep;
                i0 = it;
            }
            continue;
        }
        i3 = i2+iinc;
        if (i3>=imax) i3 = i3-imax;

        /* if factor is 4 */
        if (ifac==4) {
            if (mu==1)
                c1 = 1.0;
            else
                c1 = -1.0;
            for (l=0; l<m; l++) {
                j00 = i0;
                j01 = i1;
                j2 = i2;
                j3 = i3;
                for (jt=0; jt<nt; jt++) {
                    t1r = z[j00]+z[j2];
                    t1i = z[j00+1]+z[j2+1];
                    t2r = z[j01]+z[j3];
                    t2i = z[j01+1]+z[j3+1];
                    y1r = z[j00]-z[j2];
                    y1i = z[j00+1]-z[j2+1];
                    y3r = c1*(z[j01]-z[j3]);
                    y3i = c1*(z[j01+1]-z[j3+1]);
                    z[j00] = t1r+t2r;
                    z[j00+1] = t1i+t2i;
                    z[j01] = y1r-y3i;
                    z[j01+1] = y1i+y3r;
                    z[j2] = t1r-t2r;
                    z[j2+1] = t1i-t2i;
                    z[j3] = y1r+y3i;
                    z[j3+1] = y1i-y3r;
                    j00 += jstep;
                    j01 += jstep;
                    j2 += jstep;
                    j3 += jstep;
                }
                it = i3+istep;
                i3 = i2+istep;
                i2 = i1+istep;
                i1 = i0+istep;
                i0 = it;
            }
            continue;
        }
        i4 = i3+iinc;
        if (i4>=imax) i4 = i4-imax;

        /* if factor is 5 */