chbevl.c

/*							chbevl.c
 *
 *	Evaluate Chebyshev series
 *
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * double x, y, coef[N], chebevl();
 *
 * y = chbevl( x, coef, N );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates the series
 *
 *        N-1
 *         - '
 *  y  =   >   coef[i] T (x/2)
 *         -            i
 *        i=0
 *
 * of Chebyshev polynomials Ti at argument x/2.
 *
 * Coefficients are stored in reverse order, i.e. the zero
 * order term is last in the array.  Note N is the number of
 * coefficients, not the order.
 *
 * If coefficients are for the interval a to b, x must
 * have been transformed to x -> 2(2x - b - a)/(b-a) before
 * entering the routine.  This maps x from (a, b) to (-1, 1),
 * over which the Chebyshev polynomials are defined.
 *
 * If the coefficients are for the inverted interval, in
 * which (a, b) is mapped to (1/b, 1/a), the transformation
 * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
 * this becomes x -> 4a/x - 1.
 *
 *
 *
 * SPEED:
 *
 * Taking advantage of the recurrence properties of the
 * Chebyshev polynomials, the routine requires one more
 * addition per loop than evaluating a nested polynomial of
 * the same degree.
 *
 */
/*							chbevl.c	*/

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1985, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140

Some software in this archive may be from the book _Methods and
Programs for Mathematical Functions_ (Prentice-Hall, 1989) or
from the Cephes Mathematical Library, a commercial product. In
either event, it is copyrighted by the author.  What you see here
may be used freely but it comes with no support or guarantee.

The two known misprints in the book are repaired here in the
source listings for the gamma function and the incomplete beta
integral.


Stephen L. Moshier
moshier@world.std.com
*/

double chbevl( x, array, n )
double x;
double array[];
int n;
{
double b0, b1, b2, *p;
int i;

p = array;
b0 = *p++;
b1 = 0.0;
i = n - 1;

do
	{
	b2 = b1;
	b1 = b0;
	b0 = x * b1  -  b2  + *p++;
	}
while( --i );

return( 0.5*(b0-b2) );
}