/* fft.f -- translated by f2c (version 19950110).  */

typedef long int integer;
typedef float real;
typedef double doublereal;
typedef struct { double r, i; } doublecomplex;
typedef long int logical;

double d_sign(double *a, double *b) {
    double x;
    x = (*a >= 0 ? *a : - *a);
    return(*b >= 0 ? x : -x);
}

/* Subroutine */ int fft_(x, n, k)
doublecomplex *x;
integer *n, *k;
{
    /* Initialized data */

    static doublereal pi2 = 6.283185307179586477;
    static doublereal gain = 1.;
    static integer no = 0;
    static integer ko = 0;

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2;

    /* Builtin functions */
    double cos(), sin(), d_sign();

    /* Local variables */
    static integer i, j, l, p, q, s;
    static doublecomplex t, u[16], v, w;
    static integer stage, index, jndex, itemp, k0;
    static doublecomplex xtemp;
    static doublereal im, re;
    static integer l2n;
    static doublereal ang;
    static logical new__;
    static integer sby2;

/* 	FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X */
/* 	(A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical */
/*       methods for engineering applications. */

/* 	X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM. */
/* 	N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<32769). */
/* 	K = 1 FOR FORWARD TRANSFORM. */
/* 	K = -1 FOR INVERSE TRANSFORM. */

    /* Parameter adjustments */
    --x;

    /* Function Body */

/* 	TEST FIRST CALL? */

    new__ = no != *n;
    if (! new__) {
	goto L2;
    }

/* 	IF FIRST CALL COMPUTE LOG2 (N). */

    l2n = 0;
    no = 1;
L1:
    ++l2n;
    no += no;
    if (no < *n) {
	goto L1;
    }
    gain = (float)1. / *n;
    ang = pi2 * gain;
    re = cos(ang);
    im = sin(ang);

/* 	COMPUTE COMPLEX EXPONENTIALS IF NOT FIRST CALL */

L2:
    if (! new__ && *k * ko >= 1) {
	goto L4;
    }
    d__1 = (real) (*k);
    d__2 = -d_sign(&im, &d__1);
    z__1.r = re, z__1.i = d__2;
    u[0].r = z__1.r, u[0].i = z__1.i;
    i__1 = l2n;
    for (i = 2; i <= i__1; ++i) {
	i__2 = i - 1;
	i__3 = i - 2;
	i__4 = i - 2;
	z__1.r = u[i__3].r * u[i__4].r - u[i__3].i * u[i__4].i, z__1.i = u[
		i__3].r * u[i__4].i + u[i__3].i * u[i__4].r;
	u[i__2].r = z__1.r, u[i__2].i = z__1.i;
/* L3: */
    }
    k0 = *k;

/* 	MAIN LOOP */

L4:
    sby2 = *n;
    i__1 = l2n;
    for (stage = 1; stage <= i__1; ++stage) {
	i__2 = stage - 1;
	v.r = u[i__2].r, v.i = u[i__2].i;
	w.r = (float)1., w.i = (float)0.;
	s = sby2;
	sby2 = s / 2;
	i__2 = sby2;
	for (l = 1; l <= i__2; ++l) {
	    i__3 = *n;
	    i__4 = s;
	    for (i = 1; i__4 < 0 ? i >= i__3 : i <= i__3; i += i__4) {
		p = i + l - 1;
		q = p + sby2;
		i__5 = p;
		i__6 = q;
		z__1.r = x[i__5].r + x[i__6].r, z__1.i = x[i__5].i + x[i__6]
			.i;
		t.r = z__1.r, t.i = z__1.i;
		i__5 = q;
		i__6 = p;
		i__7 = q;
		z__2.r = x[i__6].r - x[i__7].r, z__2.i = x[i__6].i - x[i__7]
			.i;
		z__1.r = z__2.r * w.r - z__2.i * w.i, z__1.i = z__2.r * w.i + 
			z__2.i * w.r;
		x[i__5].r = z__1.r, x[i__5].i = z__1.i;
		i__5 = p;
		x[i__5].r = t.r, x[i__5].i = t.i;
/* L5: */
	    }
	    z__1.r = w.r * v.r - w.i * v.i, z__1.i = w.r * v.i + w.i * v.r;
	    w.r = z__1.r, w.i = z__1.i;
/* L6: */
	}
/* L7: */
    }

/* 	REORDER THE ELEMENTS BY BIT REVERSAL */

    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	index = i - 1;
	jndex = 0;
	i__2 = l2n;
	for (j = 1; j <= i__2; ++j) {
	    jndex += jndex;
	    itemp = index / 2;
	    if (itemp + itemp != index) {
		++jndex;
	    }
	    index = itemp;
/* L8: */
	}
	j = jndex + 1;
	if (j < i) {
	    goto L9;
	}
	i__2 = j;
	xtemp.r = x[i__2].r, xtemp.i = x[i__2].i;
	i__2 = j;
	i__4 = i;
	x[i__2].r = x[i__4].r, x[i__2].i = x[i__4].i;
	i__2 = i;
	x[i__2].r = xtemp.r, x[i__2].i = xtemp.i;
L9:
	;
    }

/* 	FORWARD TRANSFORM DONE */

    if (*k > 0) {
	return 0;
    }

/* 	INVERSE TRANSFORM */

    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	i__2 = i;
	i__4 = i;
	z__1.r = gain * x[i__4].r, z__1.i = gain * x[i__4].i;
	x[i__2].r = z__1.r, x[i__2].i = z__1.i;
/* L10: */
    }
    return 0;
} /* fft_ */