fftpack.cpp

Add c++ support for Gaussian Quadrature v1.1

/* 
 * 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 ([e-mail]);
 * 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); */

#endif

static 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;