/*  * This file is based largely on the following software distribution: *  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  *                              FFTPACK *  * Reference                                                                                                                         *    P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations *    (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.                                                                                                                    *  *     http://www.netlib.org/fftpack/ *  * Updated to single, double, and extended precision, * and translated to ISO-Standard C/C++ (without aliasing) * on 10 October 2005 by Andrew Fernandes <andrew_AT_fernandes.org> *  *                   Version 4  April 1985 *  *      A Package of Fortran Subprograms for the Fast Fourier *       Transform of Periodic and other Symmetric Sequences *  *                          by *  *                   Paul N Swarztrauber *  *   National Center for Atmospheric Research, Boulder, Colorado 80307, *  *    which is sponsored by the National Science Foundation *  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * There appears to be no explicit license for FFTPACK. However, the * package has been incorporated verbatim into a large number of software * systems over the years with numerous types of license without complaint * from the original author; therefore it would appear * that the code is effectively public domain. If you are in doubt, * however, you will need to contact the author or the  National Center * for Atmospheric Research to be sure. *  * All the changes from the original FFTPACK to the current file * fall under the following BSD-style open-source license: *  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * Copyright (c) 2005, Andrew Fernandes (andrew@fernandes.org); * All rights reserved. *   * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: *   * - Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. *  * - Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. *  * - Neither the name of the North Carolina State University nor the * names of its contributors may be used to endorse or promote products * derived from this software without specific prior written permission. *   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. * */#include "fftpack.h"#ifdef __cplusplus#include <cmath> /* the correct precision will be automatically selected */using std::cos;using std::sin;#else /* ! __cplusplus */#include <math.h> /* you must define/typedef the functions 'cos/cosf/cosl' and 'sin/sinf/sinl' as appropriate *//* real_t cos(real_t); *//* real_t sin(real_t); */#endifstatic void cfftb1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );static void cfftf1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );static void cffti1( integer_t *n, real_t *wa, integer_t *ifac );static void cosqb1( integer_t *n, real_t *x, real_t *w, real_t *xh, integer_t *ifac );static void cosqf1( integer_t *n, real_t *x, real_t *w, real_t *xh, integer_t *ifac );static void ezfft1( integer_t *n, real_t *wa, integer_t *ifac );static void passb( integer_t *nac, integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );static void passb2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );static void passb3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );static void passb4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );static void passb5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );static void passf( integer_t *nac, integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );static void passf2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );static void passf3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );static void passf4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );static void passf5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );static void radb2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );static void radb3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );static void radb4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );static void radb5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );static void radbg( integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );static void radf2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );static void radf3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );static void radf4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );static void radf5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );static void radfg( integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );static void rfftb1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );static void rfftf1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );static void rffti1( integer_t *n, real_t *wa, integer_t *ifac );static void sint1( integer_t *n, real_t *war, real_t *was, real_t *xh, real_t *x, integer_t *ifac );/* Subroutine */ void cfftb(integer_t *n, real_t *c__, real_t *wsave, 	integer_t *ifac){    integer_t iw1;    /* Parameter adjustments */    --ifac;    --wsave;    --c__;    /* Function Body */    if (*n == 1) {	return;    }    iw1 = *n + *n + 1;    cfftb1(n, &c__[1], &wsave[1], &wsave[iw1], &ifac[1]);    return;} /* cfftb_ *//* Subroutine */ static void cfftb1(integer_t *n, real_t *c__, real_t *ch, 	real_t *wa, integer_t *ifac){    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t i__, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, 	    idl1, idot;    /* Parameter adjustments */    --ifac;    --wa;    --ch;    --c__;    /* Function Body */    nf = ifac[2];    na = 0;    l1 = 1;    iw = 1;    i__1 = nf;    for (k1 = 1; k1 <= i__1; ++k1) {	ip = ifac[k1 + 2];	l2 = ip * l1;	ido = *n / l2;	idot = ido + ido;	idl1 = idot * l1;	if (ip != 4) {	    goto L103;	}	ix2 = iw + idot;	ix3 = ix2 + idot;	if (na != 0) {	    goto L101;	}	passb4(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);	goto L102;L101:	passb4(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);L102:	na = 1 - na;	goto L115;L103:	if (ip != 2) {	    goto L106;	}	if (na != 0) {	    goto L104;	}	passb2(&idot, &l1, &c__[1], &ch[1], &wa[iw]);	goto L105;L104:	passb2(&idot, &l1, &ch[1], &c__[1], &wa[iw]);L105:	na = 1 - na;	goto L115;L106:	if (ip != 3) {	    goto L109;	}	ix2 = iw + idot;	if (na != 0) {	    goto L107;	}	passb3(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);	goto L108;L107:	passb3(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);L108:	na = 1 - na;	goto L115;L109:	if (ip != 5) {	    goto L112;	}	ix2 = iw + idot;	ix3 = ix2 + idot;	ix4 = ix3 + idot;	if (na != 0) {	    goto L110;	}	passb5(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[		ix4]);	goto L111;L110:	passb5(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[		ix4]);L111:	na = 1 - na;	goto L115;L112:	if (na != 0) {	    goto L113;	}	passb(&nac, &idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1]		, &ch[1], &wa[iw]);	goto L114;L113:	passb(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1], 		&c__[1], &wa[iw]);L114:	if (nac != 0) {	    na = 1 - na;	}L115:	l1 = l2;	iw += (ip - 1) * idot;/* L116: */    }    if (na == 0) {	return;    }    n2 = *n + *n;    i__1 = n2;    for (i__ = 1; i__ <= i__1; ++i__) {	c__[i__] = ch[i__];/* L117: */    }    return;} /* cfftb1_ *//* Subroutine */ void cfftf(integer_t *n, real_t *c__, real_t *wsave, 	integer_t *ifac){    integer_t iw1;    /* Parameter adjustments */    --ifac;    --wsave;    --c__;    /* Function Body */    if (*n == 1) {	return;    }    iw1 = *n + *n + 1;    cfftf1(n, &c__[1], &wsave[1], &wsave[iw1], &ifac[1]);    return;} /* cfftf_ *//* Subroutine */ static void cfftf1(integer_t *n, real_t *c__, real_t *ch, 	real_t *wa, integer_t *ifac){    /* System generated locals */    integer_t i__1;    /* Local variables */    integer_t i__, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, 	    idl1, idot;    /* Parameter adjustments */    --ifac;    --wa;    --ch;    --c__;    /* Function Body */    nf = ifac[2];    na = 0;    l1 = 1;    iw = 1;    i__1 = nf;    for (k1 = 1; k1 <= i__1; ++k1) {	ip = ifac[k1 + 2];	l2 = ip * l1;	ido = *n / l2;	idot = ido + ido;	idl1 = idot * l1;	if (ip != 4) {	    goto L103;	}	ix2 = iw + idot;	ix3 = ix2 + idot;	if (na != 0) {	    goto L101;	}	passf4(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);	goto L102;L101:	passf4(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);L102:	na = 1 - na;	goto L115;L103:	if (ip != 2) {	    goto L106;	}	if (na != 0) {	    goto L104;	}	passf2(&idot, &l1, &c__[1], &ch[1], &wa[iw]);	goto L105;L104:	passf2(&idot, &l1, &ch[1], &c__[1], &wa[iw]);L105:	na = 1 - na;	goto L115;L106:	if (ip != 3) {	    goto L109;	}	ix2 = iw + idot;	if (na != 0) {	    goto L107;	}	passf3(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);	goto L108;L107:	passf3(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);L108:	na = 1 - na;	goto L115;L109:	if (ip != 5) {	    goto L112;	}	ix2 = iw + idot;	ix3 = ix2 + idot;	ix4 = ix3 + idot;	if (na != 0) {	    goto L110;	}	passf5(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[		ix4]);	goto L111;L110:	passf5(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[		ix4]);L111:	na = 1 - na;	goto L115;L112:	if (na != 0) {	    goto L113;	}	passf(&nac, &idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1]		, &ch[1], &wa[iw]);	goto L114;L113:	passf(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1], 		&c__[1], &wa[iw]);L114:	if (nac != 0) {	    na = 1 - na;	}L115:	l1 = l2;	iw += (ip - 1) * idot;/* L116: */    }    if (na == 0) {	return;    }    n2 = *n + *n;    i__1 = n2;    for (i__ = 1; i__ <= i__1; ++i__) {	c__[i__] = ch[i__];/* L117: */    }    return;} /* cfftf1_ *//* Subroutine */ void cffti(integer_t *n, real_t *wsave, integer_t *ifac){    integer_t iw1;    /* Parameter adjustments */    --ifac;    --wsave;    /* Function Body */    if (*n == 1) {	return;    }    iw1 = *n + *n + 1;    cffti1(n, &wsave[iw1], &ifac[1]);    return;} /* cffti_ *//* Subroutine */ static void cffti1(integer_t *n, real_t *wa, integer_t *ifac){    /* Initialized data */    static integer_t ntryh[4] = { 3,4,2,5 };    /* System generated locals */    integer_t i__1, i__2, i__3;    /* Local variables */    integer_t i__, j, i1, k1, l1, l2, ib;    real_t fi;    integer_t ld, ii, nf, ip, nl, nq, nr;    real_t arg;    integer_t ido, ipm;    real_t tpi, argh;    integer_t idot, ntry=0;    real_t argld;    /* 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;