/* gsl/poly/gsl_poly_orth.c
*
* Copyright (C) 2006 Richard J. Mathar
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include
#include
#include
/**
* @param[out] dd divided differences for Hermite Polynomial H(n,x).
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[out] xa abscissa values for the dd array
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[in] size One more than the index n of the polynomial; also the
* length of the arrays dd and xa.
* A call of the function prepares the arrays dd and xa for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthH (double dd[], double xa[], const size_t size)
{
int l ;
const int n=size-1 ;
double fact= gsl_pow_int(2.,n) ;
/* We use the representation 25.1.10 of the divided differences,
* evaluating the derivatives of the Hermite polynomial at x=0.
*/
for(l=0 ; l < size ; l++)
xa[l]=0. ;
/* Since the polynomial has even parity for even n and odd parity for
* odd n, the l'th derivative at x=0 is zero for odd n-l.
*/
for(l= n-1; l >= 0 ; l -= 2)
dd[l]=0. ;
/* the l'th derivative is generally computed by eq 22.3.10 of
* Abramowitz & Stegun
*/
for(l= n; l >= 0 ; l -= 2)
{
dd[l]= fact ;
fact *= -0.5*l*(l-1)/(double)(n-l+2) ;
}
return GSL_SUCCESS;
}
/**
* @param[out] dd divided differences for Laguerre Polynomial L(n,x).
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[out] xa abscissa values for the dd array
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[in] size One more than the index n of the polynomial; also the
* length of the arrays dd and xa.
* After a call of the function, the arrays dd and xa can be used for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthL (double dd[], double xa[], const size_t size)
{
int l ;
const int n=size-1 ;
double fact= 1.0 ;
/* We use the representation 25.1.10 of the divided differences,
* evaluating the derivatives of the Laguerre polynomial at x=0.
*/
for(l=0 ; l < size ; l++)
xa[l]=0. ;
/* the l'th derivative is generally computed by eq 22.3.9 of
* Abramowitz & Stegun for the case alpha=0.
*/
for(l= 0; l <= n ; l ++)
{
dd[l]= fact ;
fact *= (l-n)/gsl_pow_2(1+l) ;
}
return GSL_SUCCESS;
}
/**
* @param[out] dd divided differences for Chebyshev Polynomial of the 1st kind T(n,x).
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[out] xa abscissa values for the dd array
* Array must have been allocated with n+1 elements on input and is modified
* on output.
* @param[in] size One more than the index n of the polynomial; also the
* length of the arrays dd and xa.
* After a call of the function, the arrays dd and xa can be used for calls
* with gsl_poly_dd_eval() or gsl_poly_dd_taylor() .
* @since 2006-09-09
* @author Richard J. Mathar
*/
int gsl_poly_dd_init_orthT (double dd[], double xa[], const size_t size)
{
int l ;
const int n=size-1 ;
double fact= gsl_pow_int(2.,n-1) ;
/* We use the representation 25.1.10 of the divided differences,
* evaluating the derivatives of the Chebyshev polynomial at x=0.
*/
for(l=0 ; l < size ; l++)
xa[l]=0. ;
/* Since the polynomial has even parity for even n and odd parity for
* odd n, the l'th derivative at x=0 is zero for odd n-l.
*/
for(l= n-1; l >= 0 ; l -= 2)
dd[l]=0. ;
/* the l'th derivative is generally computed by eq 22.3.6 of
* Abramowitz & Stegun, here only at x=0.
*/
for(l= n; l >= 0 ; l -= 2)
{
dd[l]= fact ;
fact *= l*(1-l)/(double)((n-l+2)*(n+l-2)) ;
}
return GSL_SUCCESS;
}