pfacc.c

seismic software,very useful

/* Copyright (c) Colorado School of Mines, 1990.
/* All rights reserved.                       */

/*
FUNCTION:  prime factor fft:  complex to complex transform, in place

PARAMETERS:
isign		i sign of isign is the sign of exponent in fourier kernel
n			i length of transform (see notes below)
z			b array of complex numbers to be transformed in place
 
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:  I. D. Hale, Colorado School of Mines, 04/27/89
*/

#define P120 0.120536680
#define P142 0.142314838
#define P173 0.173648178
#define P222 0.222520934
#define P239 0.239315664
#define P281 0.281732557
#define P342 0.342020143
#define P354 0.354604887
#define P382 0.382683432
#define P415 0.415415013
#define P433 0.433883739
#define P464 0.464723172
#define P540 0.540640817
#define P559 0.559016994
#define P568 0.568064747
#define P587 0.587785252
#define P623 0.623489802
#define P642 0.642787610
#define P654 0.654860734
#define P663 0.663122658
#define P707 0.707106781
#define P748 0.748510748
#define P755 0.755749574
#define P766 0.766044443
#define P781 0.781831482
#define P822 0.822983866
#define P841 0.841253533
#define P866 0.866025404
#define P885 0.885456026
#define P900 0.900968868
#define P909 0.909631995
#define P923 0.923879533
#define P935 0.935016243
#define P939 0.939692621
#define P951 0.951056516
#define P959 0.959492974
#define P970 0.970941817
#define P974 0.974927912
#define P984 0.984807753
#define P989 0.989821442
#define P992 0.992708874
#define NFAX 10

typedef struct _complexStruct {
	float r,i;
} complex;

void pfacc (int isign, int n, complex cz[])
{
	static int kfax[] = { 16,13,11,9,8,7,5,4,3,2 };
	register float *z=(float*)cz;
	register int j0,j1,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,mu,l;
	float 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;
        j0 = 0;
        j1 = j0+jinc;

		/* if factor is 2 */
        if (ifac==2) {
			for (l=0; l<m; l++) {
				t1r = z[j0]-z[j1];
				t1i = z[j0+1]-z[j1+1];
				z[j0] = z[j0]+z[j1];
				z[j0+1] = z[j0+1]+z[j1+1];
				z[j1] = t1r;
				z[j1+1] = t1i;
				jt = j1+2;
				j1 = j0+2;
				j0 = jt;
			}
			continue;
		}
        j2 = j1+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[j1]+z[j2];
				t1i = z[j1+1]+z[j2+1];
				y1r = z[j0]-0.5*t1r;
				y1i = z[j0+1]-0.5*t1i;
				y2r = c1*(z[j1]-z[j2]);
				y2i = c1*(z[j1+1]-z[j2+1]);
				z[j0] = z[j0]+t1r;
				z[j0+1] = z[j0+1]+t1i;
				z[j1] = y1r-y2i;
				z[j1+1] = y1i+y2r;
				z[j2] = y1r+y2i;
				z[j2+1] = y1i-y2r;
				jt = j2+2;
				j2 = j1+2;
				j1 = j0+2;
				j0 = 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[j0]+z[j2];
				t1i = z[j0+1]+z[j2+1];
				t2r = z[j1]+z[j3];
				t2i = z[j1+1]+z[j3+1];
				y1r = z[j0]-z[j2];
				y1i = z[j0+1]-z[j2+1];
				y3r = c1*(z[j1]-z[j3]);
				y3i = c1*(z[j1+1]-z[j3+1]);
				z[j0] = t1r+t2r;
				z[j0+1] = t1i+t2i;
				z[j1] = y1r-y3i;
				z[j1+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 = j1+2;
				j1 = j0+2;
				j0 = 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[j1]+z[j4];
				t1i = z[j1+1]+z[j4+1];
				t2r = z[j2]+z[j3];
				t2i = z[j2+1]+z[j3+1];
				t3r = z[j1]-z[j4];
				t3i = z[j1+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[j0]-0.25*t5r;
				t7i = z[j0+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[j0] = z[j0]+t5r;
				z[j0+1] = z[j0+1]+t5i;
				z[j1] = y1r-y4i;
				z[j1+1] = y1i+y4r;
				z[j2] = y2r-y3i;
				z[j2+1] = y2i+y3r;
				z[j3] = y2r+y3i;
				z[j3+1] = y2i-y3r;
				z[j4] = y1r+y4i;
				z[j4+1] = y1i-y4r;
				jt = j4+2;
				j4 = j3+2;
				j3 = j2+2;
				j2 = j1+2;
				j1 = j0+2;
				j0 = jt;
			}
			continue;
		}
		j5 = j4+jinc;
		if (j5>=jmax) j5 = j5-jmax;
		j6 = j5+jinc;
		if (j6>=jmax) j6 = j6-jmax;

		/* if factor is 7 */
		if (ifac==7) {
			if (mu==1) {
				c1 = P623;
				c2 = -P222;
				c3 = -P900;
				c4 = P781;
				c5 = P974;
				c6 = P433;
			} else if (mu==2) {
				c1 = -P222;
				c2 = -P900;
				c3 = P623;
				c4 = P974;
				c5 = -P433;
				c6 = -P781;
			} else if (mu==3) {
				c1 = -P900;
				c2 = P623;
				c3 = -P222;
				c4 = P433;
				c5 = -P781;
				c6 = P974;
			} else if (mu==4) {
				c1 = -P900;
				c2 = P623;
				c3 = -P222;
				c4 = -P433;
				c5 = P781;
				c6 = -P974;
			} else if (mu==5) {
				c1 = -P222;
				c2 = -P900;
				c3 = P623;
				c4 = -P974;
				c5 = P433;
				c6 = P781;
			} else {
				c1 = P623;
				c2 = -P222;
				c3 = -P900;
				c4 = -P781;
				c5 = -P974;
				c6 = -P433;
			}
			for (l=0; l<m; l++) {
				t1r = z[j1]+z[j6];
				t1i = z[j1+1]+z[j6+1];
				t2r = z[j2]+z[j5];
				t2i = z[j2+1]+z[j5+1];
				t3r = z[j3]+z[j4];
				t3i = z[j3+1]+z[j4+1];
				t4r = z[j1]-z[j6];
				t4i = z[j1+1]-z[j6+1];
				t5r = z[j2]-z[j5];
				t5i = z[j2+1]-z[j5+1];
				t6r = z[j3]-z[j4];
				t6i = z[j3+1]-z[j4+1];
				t7r = z[j0]-0.5*t3r;
				t7i = z[j0+1]-0.5*t3i;
				t8r = t1r-t3r;
				t8i = t1i-t3i;
				t9r = t2r-t3r;
				t9i = t2i-t3i;
				y1r = t7r+c1*t8r+c2*t9r;
				y1i = t7i+c1*t8i+c2*t9i;
				y2r = t7r+c2*t8r+c3*t9r;
				y2i = t7i+c2*t8i+c3*t9i;
				y3r = t7r+c3*t8r+c1*t9r;
				y3i = t7i+c3*t8i+c1*t9i;
				y4r = c6*t4r-c4*t5r+c5*t6r;
				y4i = c6*t4i-c4*t5i+c5*t6i;
				y5r = c5*t4r-c6*t5r-c4*t6r;
				y5i = c5*t4i-c6*t5i-c4*t6i;
				y6r = c4*t4r+c5*t5r+c6*t6r;
				y6i = c4*t4i+c5*t5i+c6*t6i;
				z[j0] = z[j0]+t1r+t2r+t3r;
				z[j0+1] = z[j0+1]+t1i+t2i+t3i;
				z[j1] = y1r-y6i;
				z[j1+1] = y1i+y6r;
				z[j2] = y2r-y5i;
				z[j2+1] = y2i+y5r;
				z[j3] = y3r-y4i;
				z[j3+1] = y3i+y4r;
				z[j4] = y3r+y4i;
				z[j4+1] = y3i-y4r;
				z[j5] = y2r+y5i;
				z[j5+1] = y2i-y5r;
				z[j6] = y1r+y6i;
				z[j6+1] = y1i-y6r;
				jt = j6+2;
				j6 = j5+2;
				j5 = j4+2;
				j4 = j3+2;
				j3 = j2+2;
				j2 = j1+2;
				j1 = j0+2;
				j0 = jt;
			}
			continue;
		}
		j7 = j6+jinc;
		if (j7>=jmax) j7 = j7-jmax;

		/* if factor is 8 */
		if (ifac==8) {
			if (mu==1) {
				c1 = 1.0;
				c2 = P707;
			} else if (mu==3) {
				c1 = -1.0;
				c2 = -P707;
			} else if (mu==5) {
				c1 = 1.0;
				c2 = -P707;
			} else {
				c1 = -1.0;
				c2 = P707;
			}
			c3 = c1*c2;
			for (l=0; l<m; l++) {
				t1r = z[j0]+z[j4];
				t1i = z[j0+1]+z[j4+1];
				t2r = z[j0]-z[j4];
				t2i = z[j0+1]-z[j4+1];
				t3r = z[j1]+z[j5];
				t3i = z[j1+1]+z[j5+1];
				t4r = z[j1]-z[j5];
				t4i = z[j1+1]-z[j5+1];
				t5r = z[j2]+z[j6];
				t5i = z[j2+1]+z[j6+1];
				t6r = c1*(z[j2]-z[j6]);
				t6i = c1*(z[j2+1]-z[j6+1]);
				t7r = z[j3]+z[j7];
				t7i = z[j3+1]+z[j7+1];
				t8r = z[j3]-z[j7];
				t8i = z[j3+1]-z[j7+1];
				t9r = t1r+t5r;
				t9i = t1i+t5i;
				t10r = t3r+t7r;
				t10i = t3i+t7i;
				t11r = c2*(t4r-t8r);
				t11i = c2*(t4i-t8i);
				t12r = c3*(t4r+t8r);
				t12i = c3*(t4i+t8i);
				y1r = t2r+t11r;
				y1i = t2i+t11i;
				y2r = t1r-t5r;
				y2i = t1i-t5i;
				y3r = t2r-t11r;
				y3i = t2i-t11i;
				y5r = t12r-t6r;
				y5i = t12i-t6i;
				y6r = c1*(t3r-t7r);
				y6i = c1*(t3i-t7i);
				y7r = t12r+t6r;
				y7i = t12i+t6i;
				z[j0] = t9r+t10r;
				z[j0+1] = t9i+t10i;
				z[j1] = y1r-y7i;
				z[j1+1] = y1i+y7r;
				z[j2] = y2r-y6i;
				z[j2+1] = y2i+y6r;
				z[j3] = y3r-y5i;
				z[j3+1] = y3i+y5r;
				z[j4] = t9r-t10r;
				z[j4+1] = t9i-t10i;
				z[j5] = y3r+y5i;
				z[j5+1] = y3i-y5r;
				z[j6] = y2r+y6i;
				z[j6+1] = y2i-y6r;
				z[j7] = y1r+y7i;
				z[j7+1] = y1i-y7r;
				jt = j7+2;
				j7 = j6+2;
				j6 = j5+2;
				j5 = j4+2;
				j4 = j3+2;
				j3 = j2+2;
				j2 = j1+2;
				j1 = j0+2;
				j0 = jt;
			}
			continue;
		}
		j8 = j7+jinc;
		if (j8>=jmax) j8 = j8-jmax;

		/* if factor is 9 */
		if (ifac==9) {
			if (mu==1) {
				c1 = P866;
				c2 = P766;
				c3 = P642;
				c4 = P173;
				c5 = P984;
			} else if (mu==2) {
				c1 = -P866;
				c2 = P173;
				c3 = P984;
				c4 = -P939;
				c5 = P342;
			} else if (mu==4) {
				c1 = P866;
				c2 = -P939;
				c3 = P342;
				c4 = P766;
				c5 = -P642;
			} else if (mu==5) {
				c1 = -P866;
				c2 = -P939;
				c3 = -P342;
				c4 = P766;
				c5 = P642;
			} else if (mu==7) {
				c1 = P866;
				c2 = P173;
				c3 = -P984;
				c4 = -P939;
				c5 = -P342;
			} else {
				c1 = -P866;
				c2 = P766;
				c3 = -P642;
				c4 = P173;
				c5 = -P984;
			}
			c6 = c1*c2;
			c7 = c1*c3;
			c8 = c1*c4;
			c9 = c1*c5;
			for (l=0; l<m; l++) {
				t1r = z[j3]+z[j6];
				t1i = z[j3+1]+z[j6+1];
				t2r = z[j0]-0.5*t1r;
				t2i = z[j0+1]-0.5*t1i;
				t3r = c1*(z[j3]-z[j6]);
				t3i = c1*(z[j3+1]-z[j6+1]);
				t4r = z[j0]+t1r;
				t4i = z[j0+1]+t1i;
				t5r = z[j4]+z[j7];
				t5i = z[j4+1]+z[j7+1];
				t6r = z[j1]-0.5*t5r;
				t6i = z[j1+1]-0.5*t5i;
				t7r = z[j4]-z[j7];
				t7i = z[j4+1]-z[j7+1];
				t8r = z[j1]+t5r;
				t8i = z[j1+1]+t5i;
				t9r = z[j2]+z[j5];
				t9i = z[j2+1]+z[j5+1];
				t10r = z[j8]-0.5*t9r;
				t10i = z[j8+1]-0.5*t9i;
				t11r = z[j2]-z[j5];
				t11i = z[j2+1]-z[j5+1];
				t12r = z[j8]+t9r;
				t12i = z[j8+1]+t9i;
				t13r = t8r+t12r;
				t13i = t8i+t12i;
				t14r = t6r+t10r;
				t14i = t6i+t10i;
				t15r = t6r-t10r;
				t15i = t6i-t10i;
				t16r = t7r+t11r;
				t16i = t7i+t11i;
				t17r = t7r-t11r;
				t17i = t7i-t11i;
				t18r = c2*t14r-c7*t17r;
				t18i = c2*t14i-c7*t17i;
				t19r = c4*t14r+c9*t17r;
				t19i = c4*t14i+c9*t17i;
				t20r = c3*t15r+c6*t16r;
				t20i = c3*t15i+c6*t16i;
				t21r = c5*t15r-c8*t16r;
				t21i = c5*t15i-c8*t16i;
				t22r = t18r+t19r;
				t22i = t18i+t19i;
				t23r = t20r-t21r;
				t23i = t20i-t21i;
				y1r = t2r+t18r;
				y1i = t2i+t18i;
				y2r = t2r+t19r;
				y2i = t2i+t19i;
				y3r = t4r-0.5*t13r;
				y3i = t4i-0.5*t13i;
				y4r = t2r-t22r;
				y4i = t2i-t22i;
				y5r = t3r-t23r;
				y5i = t3i-t23i;
				y6r = c1*(t8r-t12r);
				y6i = c1*(t8i-t12i);
				y7r = t21r-t3r;
				y7i = t21i-t3i;
				y8r = t3r+t20r;
				y8i = t3i+t20i;
				z[j0] = t4r+t13r;
				z[j0+1] = t4i+t13i;
				z[j1] = y1r-y8i;
				z[j1+1] = y1i+y8r;
				z[j2] = y2r-y7i;
				z[j2+1] = y2i+y7r;
				z[j3] = y3r-y6i;
				z[j3+1] = y3i+y6r;
				z[j4] = y4r-y5i;
				z[j4+1] = y4i+y5r;
				z[j5] = y4r+y5i;
				z[j5+1] = y4i-y5r;
				z[j6] = y3r+y6i;
				z[j6+1] = y3i-y6r;
				z[j7] = y2r+y7i;
				z[j7+1] = y2i-y7r;
				z[j8] = y1r+y8i;
				z[j8+1] = y1i-y8r;
				jt = j8+2;
				j8 = j7+2;
				j7 = j6+2;
				j6 = j5+2;
				j5 = j4+2;
				j4 = j3+2;
				j3 = j2+2;
				j2 = j1+2;
				j1 = j0+2;
				j0 = jt;
			}
			continue;
		}
		j9 = j8+jinc;
		if (j9>=jmax) j9 = j9-jmax;
		j10 = j9+jinc;
		if (j10>=jmax) j10 = j10-jmax;

		/* if factor is 11 */
		if (ifac==11) {
			if (mu==1) {
				c1 = P841;
				c2 = P415;
				c3 = -P142;
				c4 = -P654;
				c5 = -P959;
				c6 = P540;
				c7 = P909;
				c8 = P989;
				c9 = P755;
				c10 = P281;
			} else if (mu==2) {
				c1 = P415;
				c2 = -P654;
				c3 = -P959;
				c4 = -P142;
				c5 = P841;
				c6 = P909;
				c7 = P755;
				c8 = -P281;
				c9 = -P989;
				c10 = -P540;
			} else if (mu==3) {
				c1 = -P142;
				c2 = -P959;
				c3 = P415;
				c4 = P841;
				c5 = -P654;
				c6 = P989;
				c7 = -P281;
				c8 = -P909;
				c9 = P540;
				c10 = P755;
			} else if (mu==4) {
				c1 = -P654;
				c2 = -P142;
				c3 = P841;
				c4 = -P959;
				c5 = P415;