/*
** by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
** Modified: R. Mayer to work with my benchmark routines.
*/

#define union struct /* for solaris */
/*#include <stdio.h>*/
#include <math.h>

#define double float /* for run-time specializer, jll */

#ifndef PI
#define PI  3.14159265358979323846
#endif

/*
       A Duhamel-Hollman split-radix dif fft
       Ref: Electronics Letters, Jan. 5, 1984
       Complex input and output data in arrays x and y
       Length is n.
*/

int fft(np, x, y) /* argument order changed for consistency with mayer - jll */
int np ;
double x[2]; 
double y[2]; 
{
	double *px,*py;
	int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4 ;
	double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt ;
	px = x - 1; 
	py = y - 1;
	i = 2; 
	m = 1; 
	while (i < np) { 
		i = i+i; 
		m = m+1; 
	};
	n = i; 
	if (n != np) {
		for (i = np+1; i <= n; i++)  { 
			*(px + i) = 0.0F; 
			*(py + i) = 0.0F; 
		};
		printf("\nuse %d point fft",n); 
	}

	n2 = n+n;
	for (k = 1;  k <= m-1; k++ ) {
		n2 = n2 / 2; 
		n4 = n2 / 4; 
		e = 2.0F * (float)PI / n2; 
		a = 0.0F;
		for (j = 1; j<= n4 ; j++) {
			a3 = 3.0F*a; 
			cc1 = cos(a); 
			ss1 = sin(a);
			cc3 = cos(a3); 
			ss3 = sin(a3); 
			a = j*e; 
			is = j; 
			id = 2*n2;
			while ( is < n ) {
				for (i0 = is; i0 <= n-1; /*i0 = i0 + id*/) {
					i1 = i0 + n4; 
					i2 = i1 + n4; 
					i3 = i2 + n4;
					r1 = *(px+i0) - *(px+i2);
					*(px+i0) = *(px+i0) + *(px+i2);
					r2 = *(px+i1) - *(px+i3);
					*(px+i1) = *(px+i1) + *(px+i3);
					s1 = *(py+i0) - *(py+i2);
					*(py+i0) = *(py+i0) + *(py+i2);
					s2 = *(py+i1) - *(py+i3);
					*(py+i1) = *(py+i1) + *(py+i3);
					s3 = r1 - s2; 
					r1 = r1 + s2; 
					s2 = r2 - s1; 
					r2 = r2 + s1;
					*(px+i2) = r1*cc1 - s2*ss1;
					*(py+i2) = -s2*cc1 - r1*ss1;
					*(px+i3) = s3*cc3 + r2*ss3;
					*(py+i3) = r2*cc3 - s3*ss3;
					i0 = i0 + id; /* jll */
				}
				is = 2*id - n2 + j; 
				id = 4*id;
			}
		}
	}

/*
---------------------Last stage, length=2 butterfly---------------------
*/
	is = 1; 
	id = 4;
	while ( is < n) {
		for (i0 = is; i0 <= n; /*i0 = i0 + id*/) {
			i1 = i0 + 1; 
			r1 = *(px+i0);
			*(px+i0) = r1 + *(px+i1);
			*(px+i1) = r1 - *(px+i1);
			r1 = *(py+i0);
			*(py+i0) = r1 + *(py+i1);
			*(py+i1) = r1 - *(py+i1);
			i0 = i0 + id; /* jll */
		}
		is = 2*id - 1; 
		id = 4 * id; 
	}

/*
--------------------------Bit reverse counter
*/
	j = 1; 
	n1 = n - 1;
	for (i = 1; i <= n1; i++) {
		if (i < j) {
			xt = *(px+j);
			*(px+j) = *(px+i); 
			*(px+i) = xt;
			xt = *(py+j); 
			*(py+j) = *(py+i);
			*(py+i) = xt;
		}
		k = n / 2; 
		while (k < j) { 
			j = j - k; 
			k = k / 2; 
		}
		j = j + k;
	}

	/*
  for (i = 1; i<=16; i++) printf("%d  %g   %g\n",i,*(px+i),(py+i));
  */
	return(n);

}