dpfafft.c

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

{    5720, 0.034404 },
{    6006, 0.037500 },
{    6160, 0.034091 },
{    6435, 0.040214 },
{    6552, 0.037221 },
{    6864, 0.042735 },
{    6930, 0.040214 },
{    7280, 0.042980 },
{    7920, 0.045872 },
{    8008, 0.049505 },
{    8190, 0.049834 },
{    8580, 0.055762 },
{    9009, 0.057034 },
{    9240, 0.054945 },
{    9360, 0.056818 },
{   10010, 0.066667 },
{   10296, 0.065502 },
{   10920, 0.068182 },
{   11088, 0.065217 },
{   11440, 0.075000 },
{   12012, 0.078534 },
{   12870, 0.087719 },
{   13104, 0.081081 },
{   13860, 0.084270 },
{   15015, 0.102740 },
{   16016, 0.106383 },
{   16380, 0.105634 },
{   17160, 0.119048 },
{   18018, 0.123967 },
{   18480, 0.119048 },
{   20020, 0.137615 },
{   20592, 0.140187 },
{   21840, 0.154639 },
{   24024, 0.168539 },
{   25740, 0.180723 },
{   27720, 0.180723 },
{   30030, 0.220588 },
{   32760, 0.241935 },
{   34320, 0.254237 },
{   36036, 0.254237 },
{   40040, 0.288462 },
{   45045, 0.357143 },
{   48048, 0.357143 },
{   51480, 0.384615 },
{   55440, 0.384615 },
{   60060, 0.454545 },
{   65520, 0.517241 },
{   72072, 0.576923 },
{   80080, 0.625000 },
{   90090, 0.833333 },
{  102960, 0.789474 },
{  120120, 1.153846 },
{  144144, 1.153846 },
{  180180, 1.875000 },
{  240240, 2.500000 },
{  360360, 3.750000 },
{  720720, 7.500000 },
};

int npfa_d (int nmin)
/*****************************************************************************
Return smallest valid n not less than nmin for prime factor fft.
******************************************************************************
Input:
nmin		lower bound on returned value (see notes below)

Returned:	valid n for prime factor fft
******************************************************************************
Notes:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 04/28/89
Modified:  Dave Hale, Colorado School of Mines, 08/05/91
	For efficiency, use pre-computed table of valid n and costs.
*****************************************************************************/
{
	int i;
	for (i=0; i<NTAB-1 && nctab[i].n<nmin; ++i);
	return nctab[i].n;
}

int npfao_d (int nmin, int nmax)
/*****************************************************************************
Return optimal n between nmin and nmax for prime factor fft.
******************************************************************************
Input:
nmin		lower bound on returned value (see notes below)
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for prime factor fft
******************************************************************************
Notes:
The returned n will be composed of mutually prime factors from
the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.
If nmin does not exceed 720720, then the returned n will not be 
less than nmin.  The optimal n is chosen to minimize the estimated
cost of performing the fft, while satisfying the constraint, if
possible, that n not exceed nmax.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/13/89
Modified:  Dave Hale, Colorado School of Mines, 08/05/91
	For efficiency, use pre-computed table of valid n and costs.
*****************************************************************************/
{
	int i,j;
	for (i=0; i<NTAB-1 && nctab[i].n<nmin; ++i);
	for (j=i+1; j<NTAB-1 && nctab[j].n<=nmax; ++j)
		if (nctab[j].c<nctab[i].c) i = j;
	return nctab[i].n;
}

int npfar_d (int nmin)
/*****************************************************************************
Return smallest valid n not less than nmin for real-to-real_complex or 
real_complex-to-real prime factor ffts.
******************************************************************************
Input:
nmin		lower bound on returned value

Returned:	valid n for real-to-real_complex/real_complex-to-real prime factor fft
******************************************************************************
Notes:
Current implemenations of real-to-real_complex and real_complex-to-real prime 
factor ffts require that the transform length n be even and that n/2 
be a valid length for a real_complex-to-real_complex prime factor fft.  The 
value returned by npfar satisfies these conditions.  Also, see notes 
for npfa.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/16/89
*****************************************************************************/
{
    return 2*npfa_d((nmin+1)/2);
}

int npfaro_d (int nmin, int nmax)
/*****************************************************************************
Return optimal n between nmin and nmax for real-to-real_complex or 
real_complex-to-real prime factor ffts
******************************************************************************
Input:
nmin		lower bound on returned value
nmax		desired (but not guaranteed) upper bound on returned value

Returned:	valid n for real-to-real_complex/real_complex-to-real prime factor fft
******************************************************************************
Notes:
Current implemenations of real-to-real_complex and real_complex-to-real prime 
factor ffts require that the transform length n be even and that n/2 
be a valid length for a real_complex-to-real_complex prime factor fft.  The 
value returned by npfaro satisfies these conditions.  Also, see notes 
for npfao.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/16/89
*****************************************************************************/
{
    return 2*npfao_d((nmin+1)/2,(nmax+1)/2);
}

/* code for generating the next block of constants
int
main()
{

	fprintf(stdout,"%7s %4s %.20f\n","#define","P120",sin(1.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P142",sin(1.0*PI/22.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P173",sin(1.0*PI/18.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P222",sin(1.0*PI/14.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P239",sin(1.0*PI/13.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P281",sin(1.0*PI/11.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P342",sin(1.0*PI/9.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P354",sin(3.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P382",sin(1.0*PI/8.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P415",sin(3.0*PI/22.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P433",sin(1.0*PI/7.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P464",sin(2.0*PI/13.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P540",sin(2.0*PI/11.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P559",(1.0/4.0+sin(1.0*PI/10.0)));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P568",sin(5.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P587",sin(1.0*PI/5.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P623",sin(3.0*PI/14.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P642",sin(2.0*PI/9.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P654",sin(5.0*PI/22.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P663",sin(3.0*PI/13.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P707",sin(1.0*PI/4.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P748",sin(7.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P755",sin(3.0*PI/11.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P766",sin(5.0*PI/18.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P781",sin(2.0*PI/7.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P822",sin(4.0*PI/13.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P841",sin(7.0*PI/22.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P866",sin(1.0*PI/3.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P885",sin(9.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P900",sin(5.0*PI/14.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P909",sin(4.0*PI/11.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P923",sin(3.0*PI/8.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P935",sin(5.0*PI/13.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P939",sin(7.0*PI/18.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P951",sin(2.0*PI/5.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P959",sin(9.0*PI/22.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P970",sin(11.0*PI/26.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P974",sin(3.0*PI/7.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P984",sin(4.0*PI/9.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P989",sin(5.0*PI/11.0));
	fprintf(stdout,"%7s %4s %.20f\n","#define","P992",sin(6.0*PI/13.0));

	return EXIT_SUCCESS;
}
*/


#define P120 0.12053668025532304764
#define P142 0.14231483827328514358
#define P173 0.17364817766693033119
#define P222 0.22252093395631439288
#define P239 0.23931566428755773890
#define P281 0.28173255684142967104
#define P342 0.34202014332566871291
#define P354 0.35460488704253562142
#define P382 0.38268343236508978178
#define P415 0.41541501300188637957
#define P433 0.43388373911755812040
#define P464 0.46472317204376850652
#define P540 0.54064081745559755543
#define P559 0.55901699437494745126
#define P568 0.56806474673115581187
#define P587 0.58778525229247313710
#define P623 0.62348980185873348336
#define P642 0.64278760968653925190
#define P654 0.65486073394528510061
#define P663 0.66312265824079519305
#define P707 0.70710678118654746172
#define P748 0.74851074817110108128
#define P755 0.75574957435425826890
#define P766 0.76604444311897801345
#define P781 0.78183148246802980363
#define P822 0.82298386589365635224
#define P841 0.84125353283118109449
#define P866 0.86602540378443859659
#define P885 0.88545602565320991051
#define P900 0.90096886790241914600
#define P909 0.90963199535451833011
#define P923 0.92387953251128673848
#define P935 0.93501624268541483342
#define P939 0.93969262078590831688
#define P951 0.95105651629515353118
#define P959 0.95949297361449736865
#define P970 0.97094181742605201180
#define P974 0.97492791218182361934
#define P984 0.98480775301220802032
#define P989 0.98982144188093268422
#define P992 0.99270887409805397272
#define NFAX 10

void pfacc_d (int isign, int n, real_complex cz[])
/*****************************************************************************
Prime factor fft:  real_complex to real_complex transform, in place
******************************************************************************
Input:
isign		sign of isign is the sign of exponent in fourier kernel
n		length of transform (see notes below)
z		array[n] of real_complex numbers to be transformed in place

Output:
z		array[n] of real_complex numbers 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 pfa to yield meaningful results.  this
restriction implies that n is restricted to the range
	1 <= n <= 720720 (= 5*7*9*11*13*16)
******************************************************************************
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, 04/27/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,jt;
	int nleft,jfax,ifac,jfac,jinc,jmax,ndiv,m,mm=0,mu=0,l;
	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;

	/* 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 */
        jinc = 2*mm;
		jmax = 2*n;
        j00 = 0;
        j01 = j00+jinc;

		/* if factor is 2 */
        if (ifac==2) {
			for (l=0; l<m; l++) {
				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;
				jt = j01+2;
				j01 = j00+2;
				j00 = jt;
			}
			continue;
		}
        j2 = j01+jinc;
        if (j2>=jmax) j2 = j2-jmax;

		/* if factor is 3 */
        if (ifac==3) {
			if (mu==1)
				c1 = P866;
			else
				c1 = -P866;
			for (l=0; l<m; l++) {
				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;
				jt = j2+2;
				j2 = j01+2;
				j01 = j00+2;
				j00 = jt;
			}
			continue;
		}
		j3 = j2+jinc;
		if (j3>=jmax) j3 = j3-jmax;

		/* if factor is 4 */
		if (ifac==4) {
			if (mu==1)
				c1 = 1.0;
			else
				c1 = -1.0;
			for (l=0; l<m; l++) {
				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;
				jt = j3+2;
				j3 = j2+2;
				j2 = j01+2;
				j01 = j00+2;
				j00 = jt;
			}
			continue;
		}
		j4 = j3+jinc;
		if (j4>=jmax) j4 = j4-jmax;

		/* if factor is 5 */
		if (ifac==5) {
			if (mu==1) {
				c1 = P559;
				c2 = P951;
				c3 = P587;
			} else if (mu==2) {
				c1 = -P559;
				c2 = P587;
				c3 = -P951;
			} else if (mu==3) {
				c1 = -P559;
				c2 = -P587;
				c3 = P951;
			} else { 
				c1 = P559;
				c2 = -P951;
				c3 = -P587;
			}
			for (l=0; l<m; l++) {
				t1r = z[j01]+z[j4];
				t1i = z[j01+1]+z[j4+1];
				t2r = z[j2]+z[j3];
				t2i = z[j2+1]+z[j3+1];
				t3r = z[j01]-z[j4];
				t3i = z[j01+1]-z[j4+1];
				t4r = z[j2]-z[j3];
				t4i = z[j2+1]-z[j3+1];
				t5r = t1r+t2r;
				t5i = t1i+t2i;
				t6r = c1*(t1r-t2r);
				t6i = c1*(t1i-t2i);
				t7r = z[j00]-0.25*t5r;
				t7i = z[j00+1]-0.25*t5i;
				y1r = t7r+t6r;
				y1i = t7i+t6i;
				y2r = t7r-t6r;
				y2i = t7i-t6i;
				y3r = c3*t3r-c2*t4r;
				y3i = c3*t3i-c2*t4i;
				y4r = c2*t3r+c3*t4r;
				y4i = c2*t3i+c3*t4i;
				z[j00] = z[j00]+t5r;