fejer2.cpp

Add c++ support for Gaussian Quadrature v1.1

/*
 * 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 "gaussqr.h"
#include "linalg.hpp"
#include "fftpack.h"

#include <new>
using std::nothrow;

#include <cmath>
using std::cos;

gaussqr_result fejer2_abscissae( const integer_t n , real_t *z , real_t *q )
/*
	function [z,q] = fejer2(n)

	Weights q and nodes z for the n-point Fejer type-2 quadrature rules,
	calculated via the inverse (fast) real_t discrete fourier transform.

	The knots are the same knots as the Clenshaw-Curtis knots, implying that
	by doubling the number of evalulation points, you can re-use the function
	evaluations at the old points, since the new set of knots will contain the
	old set of knots. In other words, the knots for an n-point fejer2 rule
	will be contained in the (2*n+1)-point fejer2 rule.
 
	Although capable of computing knots and weights for arbitrary 'n', the
	computation will be much faster and efficient if (n+1) has factors only
	from the set from the integer_t set {2,3,4,5}.
	
	Therefore one recommended sequence of n is {3,7,15,31,63,127,255,511,1023,...}.
	Three (or five) times this set would also be acceptable.
  
*/
{
	if ( n < 3 || z == 0 || q == 0 )
		return gaussqr_illegal_argument;
	
	gaussqr_result rv = gaussqr_success;

	integer_t n_fft = n + 1;
	const real_t pi = 3.1415926535897932384626433832795028841971693993751058209749445923;

	// allocate memory
	real_t *v = new(nothrow) real_t[n_fft];
	real_t *work = new(nothrow) real_t[2*n_fft];
	integer_t *ifac = new(nothrow) integer_t[sizeof(integer_t)*8];
	if ( v == 0 || work == 0 || ifac == 0 ) {
		rv = gaussqr_memory_allocation_error;
		goto done;
	}
		
	// calculate the knots
	for ( integer_t k = 1; k < n_fft; k++ ) {
		z[k-1] = -cos((k*pi)/n_fft);
	}
	
	// calculate the transformed weights
	v[0] = 2.0;
	for ( integer_t k = 1; k < n_fft; k++ ) v[k] = 0.0;
	for ( integer_t k = 1; k < n_fft/2; k++ ) {
		v[2*k-1] = 2.0/(1.0-4.0*(k*k));
	}
	v[n_fft-1-(n_fft&1)] = (n_fft-3.0)/(2.0*(n_fft/2)-1.0)-1.0;
	
    // take the fourier transform, and scale
	rffti(&n_fft,work,ifac);
	rfftb(&n_fft,v,work,ifac);
	for ( integer_t i = 0; i < n_fft; i++ ) {
		v[i] /= n_fft;
	}
	
	// transfer u to q, removing first (zero) element
	for ( integer_t i = 1; i < n_fft; i++ ) {
		q[i-1] = v[i];
	}

done:
	delete[](ifac);
    delete[](work);
	delete[](v);
	return(rv);
}


gaussqr_result map_fejer2_domain( const real_t a, const real_t b , const domain_type type , const integer_t n , const real_t *x , real_t *y , real_t *dy )
/*
	Given an interval [a,b] and a domain_type, this subroutine maps the interval [-1,1] to [a,b],
	an interval that may be closed, left-infinite, right-infinite, or totally infinite.
 
	The input points are x, and should be enclosed in the open interval (-1:1). This restriction is NOT checked for.
 
	On output, the points y are the mapped points to the new interval, and the dy array gives the derivative of the transformation.
 
	Note that the change of integration variables is $\int_{-1}^{+1} f(x) \, dx = \int_{a}^{b} f(y(x)) dy \, dx$, and this is
	used in the quadrature approximations to integration.
 
*/
{
	if ( a >= b || n < 3 || x == 0 || y == 0 || dy == 0 )
		return(gaussqr_illegal_argument);
	
	switch(type) {
		
		case domain_finite :
			for ( integer_t i = 0; i < n; i++ ) {
				const real_t bpa = b + a;
				const real_t bma = b - a;
				y[i] = 0.5 * ( x[i]*bma + bpa );
				dy[i] = 0.5 * bma;
			}
			break;
	
		case domain_left_infinite :
			for ( integer_t i = 0; i < n; i++ ) {
				const real_t x_m1 = (1.0-x[i]);
				const real_t x_p1 = (1.0+x[i]);
				y[i] = x_m1/x_p1 + b;
				dy[i] = 2.0/(x_p1*x_p1);
			}
			break;
			
		case domain_right_infinite :
			for ( integer_t i = 0; i < n; i++ ) {
				const real_t x_m1 = (1.0-x[i]);
				const real_t x_p1 = (1.0+x[i]);
				y[i] = x_p1/x_m1 + a;
				dy[i] = 2.0/(x_m1*x_m1);
			}
			break;
			
		case domain_infinite :
			for ( integer_t i = 0; i < n; i++ ) {
				const real_t x_m1 = (1.0-x[i]);
				const real_t x_p1 = (1.0+x[i]);
				y[i] = (2.0*x[i])/(x_m1*x_p1);
				dy[i] = (2.0*(1.0+(x[i]*x[i])))/(x_m1*x_m1*x_p1*x_p1);
			}
			break;
			
		default :
			return(gaussqr_illegal_argument);
		
	}
	
	return(gaussqr_success);
}