}    i__1 = nf;    for (i__ = 2; i__ <= i__1; ++i__) {	ib = nf - i__ + 2;	ifac[ib + 2] = ifac[ib + 1];/* L106: */    }    ifac[3] = 2;L107:    if (nl != 1) {	goto L104;    }    ifac[1] = *n;    ifac[2] = nf;    tpi = REAL_CONSTANT(6.283185307179586476925286766559005768394338798750211619498891846);    argh = tpi / (real_t) (*n);    i__ = 2;    l1 = 1;    i__1 = nf;    for (k1 = 1; k1 <= i__1; ++k1) {	ip = ifac[k1 + 2];	ld = 0;	l2 = l1 * ip;	ido = *n / l2;	idot = ido + ido + 2;	ipm = ip - 1;	i__2 = ipm;	for (j = 1; j <= i__2; ++j) {	    i1 = i__;	    wa[i__ - 1] = REAL_CONSTANT(1.0);	    wa[i__] = REAL_CONSTANT(0.0);	    ld += l1;	    fi = REAL_CONSTANT(0.0);	    argld = (real_t) ld * argh;	    i__3 = idot;	    for (ii = 4; ii <= i__3; ii += 2) {		i__ += 2;		fi += REAL_CONSTANT(1.0);		arg = fi * argld;		wa[i__ - 1] = cos(arg);		wa[i__] = sin(arg);/* L108: */	    }	    if (ip <= 5) {		goto L109;	    }	    wa[i1 - 1] = wa[i__ - 1];	    wa[i1] = wa[i__];L109:	    ;	}	l1 = l2;/* L110: */    }    return;} /* cffti1_ *//* Subroutine */ void cosqb(integer_t *n, real_t *x, real_t *wsave, 	integer_t *ifac){    /* Initialized data */    static real_t tsqrt2 = 	    REAL_CONSTANT(2.82842712474619009760337744841939615713934375053896146353359476);    /* System generated locals */    integer_t i__1;    /* Local variables */    real_t x1;    /* Parameter adjustments */    --ifac;    --wsave;    --x;    /* Function Body */    if ((i__1 = *n - 2) < 0) {	goto L101;    } else if (i__1 == 0) {	goto L102;    } else {	goto L103;    }L101:    x[1] *= REAL_CONSTANT(4.0);    return;L102:    x1 = (x[1] + x[2]) * REAL_CONSTANT(4.0);    x[2] = tsqrt2 * (x[1] - x[2]);    x[1] = x1;    return;L103:    cosqb1(n, &x[1], &wsave[1], &wsave[*n + 1], &ifac[1]);    return;} /* cosqb_ *//* Subroutine */ static void cosqb1(integer_t *n, real_t *x, real_t *w, 	real_t *xh, integer_t *ifac){    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t i__, k, kc, np2, ns2;    real_t xim1;    integer_t modn;    /* Parameter adjustments */    --ifac;    --xh;    --w;    --x;    /* Function Body */    ns2 = (*n + 1) / 2;    np2 = *n + 2;    i__1 = *n;    for (i__ = 3; i__ <= i__1; i__ += 2) {	xim1 = x[i__ - 1] + x[i__];	x[i__] -= x[i__ - 1];	x[i__ - 1] = xim1;/* L101: */    }    x[1] += x[1];    modn = *n % 2;    if (modn == 0) {	x[*n] += x[*n];    }    rfftb(n, &x[1], &xh[1], &ifac[1]);    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np2 - k;	xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k];	xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc];/* L102: */    }    if (modn == 0) {	x[ns2 + 1] = w[ns2] * (x[ns2 + 1] + x[ns2 + 1]);    }    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np2 - k;	x[k] = xh[k] + xh[kc];	x[kc] = xh[k] - xh[kc];/* L103: */    }    x[1] += x[1];    return;} /* cosqb1_ *//* Subroutine */ void cosqf(integer_t *n, real_t *x, real_t *wsave, 	integer_t *ifac){    /* Initialized data */    static real_t sqrt2 = 	    REAL_CONSTANT(1.41421356237309504880168872420969807856967187536948073176679738);    /* System generated locals */    integer_t i__1;    /* Local variables */    real_t tsqx;    /* Parameter adjustments */    --ifac;    --wsave;    --x;    /* Function Body */    if ((i__1 = *n - 2) < 0) {	goto L102;    } else if (i__1 == 0) {	goto L101;    } else {	goto L103;    }L101:    tsqx = sqrt2 * x[2];    x[2] = x[1] - tsqx;    x[1] += tsqx;L102:    return;L103:    cosqf1(n, &x[1], &wsave[1], &wsave[*n + 1], &ifac[1]);    return;} /* cosqf_ *//* Subroutine */ static void cosqf1(integer_t *n, real_t *x, real_t *w, 	real_t *xh, integer_t *ifac){    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t i__, k, kc, np2, ns2;    real_t xim1;    integer_t modn;    /* Parameter adjustments */    --ifac;    --xh;    --w;    --x;    /* Function Body */    ns2 = (*n + 1) / 2;    np2 = *n + 2;    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np2 - k;	xh[k] = x[k] + x[kc];	xh[kc] = x[k] - x[kc];/* L101: */    }    modn = *n % 2;    if (modn == 0) {	xh[ns2 + 1] = x[ns2 + 1] + x[ns2 + 1];    }    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np2 - k;	x[k] = w[k - 1] * xh[kc] + w[kc - 1] * xh[k];	x[kc] = w[k - 1] * xh[k] - w[kc - 1] * xh[kc];/* L102: */    }    if (modn == 0) {	x[ns2 + 1] = w[ns2] * xh[ns2 + 1];    }    rfftf(n, &x[1], &xh[1], &ifac[1]);    i__1 = *n;    for (i__ = 3; i__ <= i__1; i__ += 2) {	xim1 = x[i__ - 1] - x[i__];	x[i__] = x[i__ - 1] + x[i__];	x[i__ - 1] = xim1;/* L103: */    }    return;} /* cosqf1_ *//* Subroutine */ void cosqi(integer_t *n, real_t *wsave, integer_t *ifac){    /* Initialized data */    static real_t pih = 	    REAL_CONSTANT(1.570796326794896619231321691639751442098584699687529104874722962);    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t k;    real_t fk, dt;    /* Parameter adjustments */    --ifac;    --wsave;    /* Function Body */    dt = pih / (real_t) (*n);    fk = REAL_CONSTANT(0.0);    i__1 = *n;    for (k = 1; k <= i__1; ++k) {	fk += REAL_CONSTANT(1.0);	wsave[k] = cos(fk * dt);/* L101: */    }    rffti(n, &wsave[*n + 1], &ifac[1]);    return;} /* cosqi_ *//* Subroutine */ void cost(integer_t *n, real_t *x, real_t *wsave, 	integer_t *ifac){    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t i__, k;    real_t c1, t1, t2;    integer_t kc;    real_t xi;    integer_t nm1, np1;    real_t x1h;    integer_t ns2;    real_t tx2, x1p3, xim2;    integer_t modn;    /* Parameter adjustments */    --ifac;    --wsave;    --x;    /* Function Body */    nm1 = *n - 1;    np1 = *n + 1;    ns2 = *n / 2;    if ((i__1 = *n - 2) < 0) {	goto L106;    } else if (i__1 == 0) {	goto L101;    } else {	goto L102;    }L101:    x1h = x[1] + x[2];    x[2] = x[1] - x[2];    x[1] = x1h;    return;L102:    if (*n > 3) {	goto L103;    }    x1p3 = x[1] + x[3];    tx2 = x[2] + x[2];    x[2] = x[1] - x[3];    x[1] = x1p3 + tx2;    x[3] = x1p3 - tx2;    return;L103:    c1 = x[1] - x[*n];    x[1] += x[*n];    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np1 - k;	t1 = x[k] + x[kc];	t2 = x[k] - x[kc];	c1 += wsave[kc] * t2;	t2 = wsave[k] * t2;	x[k] = t1 - t2;	x[kc] = t1 + t2;/* L104: */    }    modn = *n % 2;    if (modn != 0) {	x[ns2 + 1] += x[ns2 + 1];    }    rfftf(&nm1, &x[1], &wsave[*n + 1], &ifac[1]);    xim2 = x[2];    x[2] = c1;    i__1 = *n;    for (i__ = 4; i__ <= i__1; i__ += 2) {	xi = x[i__];	x[i__] = x[i__ - 2] - x[i__ - 1];	x[i__ - 1] = xim2;	xim2 = xi;/* L105: */    }    if (modn != 0) {	x[*n] = xim2;    }L106:    return;} /* cost_ *//* Subroutine */ void costi(integer_t *n, real_t *wsave, integer_t *ifac){    /* Initialized data */    static real_t pi = 	    REAL_CONSTANT(3.141592653589793238462643383279502884197169399375158209749445923);    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t k, kc;    real_t fk, dt;    integer_t nm1, np1, ns2;    /* Parameter adjustments */    --ifac;    --wsave;    /* Function Body */    if (*n <= 3) {	return;    }    nm1 = *n - 1;    np1 = *n + 1;    ns2 = *n / 2;    dt = pi / (real_t) nm1;    fk = REAL_CONSTANT(0.0);    i__1 = ns2;    for (k = 2; k <= i__1; ++k) {	kc = np1 - k;	fk += REAL_CONSTANT(1.0);	wsave[k] = sin(fk * dt) * REAL_CONSTANT(2.0);	wsave[kc] = cos(fk * dt) * REAL_CONSTANT(2.0);/* L101: */    }    rffti(&nm1, &wsave[*n + 1], &ifac[1]);    return;} /* costi_ *//* Subroutine */ static void ezfft1(integer_t *n, real_t *wa, integer_t *ifac){    /* Initialized data */    static integer_t ntryh[4] = { 4,2,3,5 };    static real_t tpi = 	    REAL_CONSTANT(6.283185307179586476925286766559005768394338798750211419498891846);    /* System generated locals */    integer_t i__1, i__2, i__3;    /* Local variables */    integer_t i__, j, k1, l1, l2, ib, ii, nf, ip, nl, is, nq, nr;    real_t ch1, sh1;    integer_t ido, ipm;    real_t dch1, ch1h, arg1, dsh1;    integer_t nfm1;    real_t argh;    integer_t ntry=0;    /* Parameter adjustments */    --ifac;    --wa;    /* Function Body */    nl = *n;    nf = 0;    j = 0;L101:    ++j;    if (j - 4 <= 0) {	goto L102;    } else {	goto L103;    }L102:    ntry = ntryh[j - 1];    goto L104;L103:    ntry += 2;L104:    nq = nl / ntry;    nr = nl - ntry * nq;    if (nr != 0) {	goto L101;    } else {	goto L105;    }L105:    ++nf;    ifac[nf + 2] = ntry;    nl = nq;    if (ntry != 2) {	goto L107;    }    if (nf == 1) {	goto L107;    }    i__1 = nf;    for (i__ = 2; i__ <= i__1; ++i__) {	ib = nf - i__ + 2;	ifac[ib + 2] = ifac[ib + 1];/* L106: */    }    ifac[3] = 2;L107:    if (nl != 1) {	goto L104;    }    ifac[1] = *n;    ifac[2] = nf;    argh = tpi / (real_t) (*n);    is = 0;    nfm1 = nf - 1;    l1 = 1;    if (nfm1 == 0) {	return;    }    i__1 = nfm1;    for (k1 = 1; k1 <= i__1; ++k1) {	ip = ifac[k1 + 2];	l2 = l1 * ip;	ido = *n / l2;