MathFunc.hpp File Reference

#include "gstlearn_export.hpp"
#include "Matrix/MatrixRectangular.hpp"

Functions
GSTLEARN_EXPORT int	mvndst_infin (double low, double sup)
GSTLEARN_EXPORT void	mvndst (int n, double *lower, double *upper, int *infin, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
GSTLEARN_EXPORT void	mvndst2n (const double *lower, const double *upper, const double *means, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
GSTLEARN_EXPORT void	mvndst4 (double *lower, double *upper, const double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
GSTLEARN_EXPORT int	besselj_table (double x, double alpha, int nb, double *b)
GSTLEARN_EXPORT double	besselj (double x, int n)
GSTLEARN_EXPORT int	besselk (double x, double alpha, int nb, double *bk)
GSTLEARN_EXPORT double	loggamma (double parameter)
GSTLEARN_EXPORT double	ut_legendre (int n, double v, bool flagNorm=true)
GSTLEARN_EXPORT VectorDouble	ut_legendreVec (int n, const VectorDouble &vecin, bool flagNorm)
GSTLEARN_EXPORT MatrixRectangular	ut_legendreMatNorm (int n, const VectorDouble &v)
GSTLEARN_EXPORT MatrixRectangular	ut_legendreAssociatedMat (int l, const VectorDouble &v, bool flagNorm=true)
GSTLEARN_EXPORT double	ut_flegendre (int n, int k0, double theta, bool flagNorm=true)
GSTLEARN_EXPORT double	ut_sphericalHarmonic (int n, int k, double theta, double phi)
GSTLEARN_EXPORT VectorDouble	ut_sphericalHarmonicVec (int n, int k, VectorDouble theta, VectorDouble phi)
GSTLEARN_EXPORT double	golden_search (double(*func_evaluate)(double test, void *user_data), void *user_data, double tolstop, double a0, double c0, double *test_loc, double *niter)
GSTLEARN_EXPORT int	ut_chebychev_count (double(*func)(double, double, const VectorDouble &), Cheb_Elem *cheb_elem, double x, const VectorDouble &blin)
GSTLEARN_EXPORT int	ut_chebychev_coeffs (double(*func)(double, double, const VectorDouble &), Cheb_Elem *cheb_elem, const VectorDouble &blin)
GSTLEARN_EXPORT void	ut_vandercorput (int n, int flag_sym, int flag_rot, int *ntri_arg, double **coor_arg)
GSTLEARN_EXPORT int	ut_icosphere (int n, int flag_rot, int *ntri_arg, double **coor_arg)
GSTLEARN_EXPORT double	ut_factorial (int k)
GSTLEARN_EXPORT void	ut_log_factorial (int nbpoly, double *factor)
GSTLEARN_EXPORT MatrixRectangular *	vanDerCorput (int n, int nd)
GSTLEARN_EXPORT MatrixRectangular	fillLegendreMatrix (const VectorDouble &r, int legendreOrder)
GSTLEARN_EXPORT int	solve_P2 (double a, double b, double c, VectorDouble &x)
GSTLEARN_EXPORT int	solve_P3 (double a, double b, double c, double d, VectorDouble &x)

Function Documentation

besselj()

GSTLEARN_EXPORT double besselj	(	double	x,
		int	n
	)

besselj_table()

GSTLEARN_EXPORT int besselj_table	(	double	x,
		double	alpha,
		int	nb,
		double *	b
	)

This routine calculates Bessel functions J SUB(NB+ALPHA) (X) for non-negative argument X, and non-negative order NB+ALPHA.

Returns

Error return code : NCALC

NCALC < -1: An argument is out of range.

1 < NCALC < NB: Not all requested function values could be calculated accurately.

BY(I) contains correct function values

Parameters

[in]	x	Working precision non-negative real argument for which J's are to be calculated.
[in]	alpha	Working precision fractional part of order for which J's are to be calculated. 0 <= ALPHA < 1.0.
[in]	nb	Integer number of functions to be calculated, NB > 0 The first function calculated is of order ALPHA, and the last is of order (NB - 1 + ALPHA).
[out]	b	Working precision output vector of length NB. If the routine terminates normally (NCALC=NB), the vector b[] contains the functions Y(ALPHA,X), ... ,Y(NB-1+ALPHA,X),

Remarks

This program is based on a program written by David J. Sookne (2)

that computes values of the Bessel functions J or I of real

argument and integer order. Modifications include the

restriction of the computation to the J Bessel function of

non-negative real argument, the extension of the computation to

arbitrary positive order, and the elimination of most underflow.

References: "A Note on Backward Recurrence Algorithms," Olver, F.

W. J., and Sookne, D. J., Math. Comp. 26, 1972, pp 941-947.

"Bessel Functions of Real Argument and Integer Order," Sookne, D.

J., NBS Jour. of Res. B. 77B, 1973, pp 125-132.

besselk()

GSTLEARN_EXPORT int besselk	(	double	x,
		double	alpha,
		int	nb,
		double *	bk
	)

This routine calculates modified Bessel functions of the second kind, K SUB(N+ALPHA) (X), for non-negative argument X and non-negative order N+ALPHA

Returns

Error return code :

NCALC < -1: An argument is out of range.

0 < NCALC < NB: Not all requested function values could be calculated accurately.

BK(I) contains correct function values for I <= NCALC, and contains the ratios

K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.

Parameters

[in]	x	Working precision non-negative real argument for which K's are to calculated. If K's are to be calculated, X must not be greater than XMAX.
[in]	alpha	Working precision fractional part of order for which K's are to be calculated. 0 <= ALPHA <1.0.
[in]	nb	Integer number of functions to be calculated, NB > 0. The first function calculated is of order ALPHA, and the last is of order (NB - 1 + ALPHA).
[out]	bk	Working precision output vector of length NB. If the routine terminates normally (NCALC=NB), the vector BK contains the functions : K(ALPHA,X), ... , K(NB-1+ALPHA,X),

Remarks

This program is based on a program written by J. B. Campbell (2)

that computes values of the Bessel functions K of real argument

and real order. References: "On Temme's Algorithm for the \remark Modified Bessel Functions of the Third Kind," Campbell, J. B.,

TOMS 6(4), Dec. 1980, pp. 581-586. "A FORTRAN IV Subroutine for

the Modified Bessel Functions of the Third Kind of Real Order and

Real Argument," Campbell, J. B., Report NRC/ERB-925, National

Research Council, Canada.

fillLegendreMatrix()

GSTLEARN_EXPORT MatrixRectangular fillLegendreMatrix	(	const VectorDouble &	r,
		int	legendreOrder
	)

golden_search()

GSTLEARN_EXPORT double golden_search	(	double(*)(double test, void *user_data)	func_evaluate,
		void *	user_data,
		double	tolstop,
		double	a0,
		double	c0,
		double *	test_loc,
		double *	niter
	)

Golden Search algorithm

Returns

The estimated value

Parameters

[in]	func_evaluate	Evaluating function
[in]	user_data	User Data
[in]	tolstop	Tolerance parameter
[in]	a0	Initial value for lower bound of interval
[in]	c0	Initial value for upper bound of interval
[out]	test_loc	Final value of the evaluating function
[out]	niter	Number of iterations

loggamma()

GSTLEARN_EXPORT double loggamma

(

double

parameter

)

Calculation of the logarithm of the gamma function

Returns

Logarithm of the gamma function

Parameters

[in]

parameter

raw value

mvndst()

GSTLEARN_EXPORT void mvndst	(	int	n,
		double *	lower,
		double *	upper,
		int *	infin,
		double *	correl,
		int	maxpts,
		double	abseps,
		double	releps,
		double *	error,
		double *	value,
		int *	inform
	)

Multivariate Normal Probability

Parameters

[in]	n	Number of variables
[in]	lower	Array of lower integration limits
[in]	upper	Array of upper integration limits
[in]	infin	Array of integration limit flags <0 for ]-Infinity; +Infinity[ =0 for ]-Infinity; upper] =1 for [lower; +Infinity[ =2 for [lower; upper]
[in]	correl	Array of correlation coefficients
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned code

Remarks

The array correl must be entered as the non-diagonal upper part

of the correlation matrix, entered by line

This subroutine uses an algorithm given in the paper:

"Numerical Computation of Multivariate Normal Probabilities", in

J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by

Alan Genz

Department of Mathematics

Washington State University

Pullman, WA 99164-3113

Email : AlanG.nosp@m.enz@.nosp@m.wsu.e.nosp@m.du

mvndst2n()

GSTLEARN_EXPORT void mvndst2n	(	const double *	lower,
		const double *	upper,
		const double *	means,
		double *	correl,
		int	maxpts,
		double	abseps,
		double	releps,
		double *	error,
		double *	value,
		int *	inform
	)

Calculate the multigaussian integral (non-normalized)

Parameters

[in]	lower	Array of lower bounds (Dimension: nvar)
[in]	upper	Array of upper bounds (Dimension: nvar)
[in]	means	Array of means (Dimension: 2)
[in]	correl	Correlation matrix (Dimension: 2*2)
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned code

mvndst4()

GSTLEARN_EXPORT void mvndst4	(	double *	lower,
		double *	upper,
		const double *	correl,
		int	maxpts,
		double	abseps,
		double	releps,
		double *	error,
		double *	value,
		int *	inform
	)

Calculate the quadri-variable gaussian integral

Parameters

[in]	lower	Array of lower bounds
[in]	upper	Array of upper bounds
[in]	correl	Correlation matrix (Dimension: 4*4)
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned code

mvndst_infin()

GSTLEARN_EXPORT int mvndst_infin	(	double	low,
		double	sup
	)

Set the flags for the bound of numerical integration

Returns

Flag for the bound

Parameters

[in]	low	Lower integration bound
[in]	sup	Upper integration bound

solve_P2()

GSTLEARN_EXPORT int solve_P2	(	double	a,
		double	b,
		double	c,
		VectorDouble &	x
	)

Find the roots of a polynomial of order 2: ax^2 + bx + c = 0

Returns

Number of real solutions

Parameters

[in]	a,b,c	Coefficients of the polynomial
[out]	x	Array of real solutions (Dimension: 2)

Remarks

When the solution is double, the returned number os 1.

solve_P3()

GSTLEARN_EXPORT int solve_P3	(	double	a,
		double	b,
		double	c,
		double	d,
		VectorDouble &	x
	)

Find the roots of a polynomial of order 3: a*x^3 + b*x^2 + c*x + d = 0

Returns

Number of real solutions

Parameters

[in]	a,b,c,d	Coefficients of the polynomial
[out]	x	Array of real solutions (Dimension: 3)

Remarks

When the solution is double, the returned number os 1.

ut_chebychev_coeffs()

GSTLEARN_EXPORT int ut_chebychev_coeffs	(	double(*)(double, double, const VectorDouble &)	func,
		Cheb_Elem *	cheb_elem,
		const VectorDouble &	blin
	)

Calculates the coefficients of the Chebychev polynomial which is an approximation of a given function

Returns

Error return code

Parameters

[in]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	blin	Array of coefficients for polynomial expansion

ut_chebychev_count()

GSTLEARN_EXPORT int ut_chebychev_count	(	double(*)(double, double, const VectorDouble &)	func,
		Cheb_Elem *	cheb_elem,
		double	x,
		const VectorDouble &	blin
	)

Evaluate the number of coefficients necessary to evaluate a function (at a sample location) at a given approximation

Returns

Minimum number of necessary coefficients

Parameters

[in]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	x	Sampling value
[in]	blin	Array of coefficients for polynomial expansion

ut_factorial()

GSTLEARN_EXPORT double ut_factorial

(

int

k

)

Calculates the factorial coefficient

Returns

Returned value

Parameters

[in]

k

Value

ut_flegendre()

GSTLEARN_EXPORT double ut_flegendre	(	int	n,
		int	k0,
		double	theta,
		bool	flagNorm
	)

Returns the Spherical Legendre normalized function. According to boost Library documentation, this returns

Y_n^m(theta, phi) = sqrt{{(2n+1)(n-m)!} / {4pi(n+m)!}} P_n^m(cos(theta))e{imphi}

with phi=0 and m = |k0|

If flagNorm is switched ON, the previous result is multiplied by:

sqrt(2 * pi)

Parameters

[in]	n	Degree
[in]	k0	Order (ABS(k0) <= n)
[in]	theta	Theta angle in radian
[in]	flagNorm	for normalized and 0 otherwise

ut_icosphere()

GSTLEARN_EXPORT int ut_icosphere	(	int	n,
		int	flag_rot,
		int *	ntri_arg,
		double **	coor_arg
	)

Generate regular Icosahedron discretization

Returns

Error return code

Parameters

[in]	n	Number of discretization steps
[in]	flag_rot	Perform a random rotation
[out]	ntri_arg	Number of points
[out]	coor_arg	Array of point coordinates (Dimension: 3*ntri)

Remarks

As random number are used in this function, a specific seed

is fixed here

ut_legendre()

GSTLEARN_EXPORT double ut_legendre	(	int	n,
		double	v,
		bool	flagNorm
	)

Returns the Legendre Function described as follows:

legendre_p(n,v) = 1/{2^n n!) d^n/dx^n [x^2-1)^n with |x|<=1

If flagNorm is switched ON, the Boost library cannot be used anymore. We have to rely on the ode provided by X; Freulon.

Parameters

[in]	n	Degree of the Legendre Polynomial to be computed (n >= 0)
[in]	v	Value for which the polynomial is evaluated (-1 <= v <= 1)
[in]	flagNorm	True for normalized and 0 otherwise

Returns

Value of the Legendre polynomial.

ut_legendreAssociatedMat()

GSTLEARN_EXPORT MatrixRectangular ut_legendreAssociatedMat	(	int	l,
		const VectorDouble &	v,
		bool	flagNorm
	)

Returns the Associated Legendre Function described as follows:

In the case flagNorm=true

Using the relations: P_0^0 (x) = 1 P_{l+1}^{l+1} (x) = - sqrt((2l+3)/(2l+2)) * (1-x^2)^{1/2} * P_{l}^{l}(x) P_{l+1}^{m} (x) = sqrt((2l+3)(2l+1)/((l-m+1)(l+m+1))) * x * P_{l}^{m}(x) - sqrt((2l+3)/(2l-1)*(l+m)/(l+m+1)*(l-m)/(l-m+1)) P_{l-1}^m(x) Changing l to l-1, P_{l}^{l} (x) = - a0 * (1-x^2)^{1/2} * P_{l-1}^{l-1}(x) with a0 = sqrt((2l+1)/(2l)) P_{l}^{m} (x) = a * x * P_{l-1}^{m}(x) - b P_{l-2}^m(x) with a = sqrt((2l+1)(2l-1)/(l-m)/(l+m)) and b = sqrt((2l+1)/(2l-3)*(l+m-1)/(l+m)*(l-m-1)/(l-m))

In the case flagNorm=false

Using the relations: P_0^0 (x) = 1 P_{l+1}^{l+1} (x) = - sqrt((2l+3)/(2l+2)) * (1-x^2)^{1/2} * P_{l}^{l}(x) P_{l+1}^{m} (x) = sqrt((2l+3)(2l+1)/((l-m+1)(l+m+1))) * x * P_{l}^{m}(x) - sqrt((2l+3)/(2l-1)*(l+m)/(l+m+1)*(l-m)/(l-m+1)) P_{l-1}^m(x) Changing l to l-1, P_{l}^{l} (x) = - a0 * (1-x^2)^{1/2} * P_{l-1}^{l-1}(x) with a0 = sqrt((2l+1)/(2l)) P_{l}^{m} (x) = a * x * P_{l-1}^{m}(x) - b P_{l-2}^m(x) with a = sqrt((2l+1)(2l-1)/(l-m)/(l+m)) and b = sqrt((2l+1)/(2l-3)*(l+m-1)/(l+m)*(l-m-1)/(l-m))

Parameters

[in]	l	Degree of the Legendre Polynomial to be computed (n >= 0)
[in]	v	Value for which the polynomial is evaluated (-1 <= v <= 1)
[in]	flagNorm	Normalized when True

Returns

A matrix of the (normalized) associated Legendre functions: P_l^m (x) for 0 <= m <= l

ut_legendreMatNorm()

GSTLEARN_EXPORT MatrixRectangular ut_legendreMatNorm	(	int	n,
		const VectorDouble &	v
	)

ut_legendreVec()

GSTLEARN_EXPORT VectorDouble ut_legendreVec	(	int	n,
		const VectorDouble &	vecin,
		bool	flagNorm
	)

ut_log_factorial()

GSTLEARN_EXPORT void ut_log_factorial	(	int	nbpoly,
		double *	factor
	)

Calculates the nbpoly log-factorial coefficients

Parameters

[in]	nbpoly	Number of terms
[out]	factor	logarithm of factorials

ut_sphericalHarmonic()

GSTLEARN_EXPORT double ut_sphericalHarmonic	(	int	n,
		int	k,
		double	theta,
		double	phi
	)

Returns the Spherical harmonic

Parameters

[in]	n	Degree of HS (n >= 0)
[in]	k	Order of the HS (-n <= k <= n)
[in]	theta	Colatitude angle in radian (0 <= theta <= pi
[in]	phi	Longitude angle in radian (0 <= phi <= 2* pi)

ut_sphericalHarmonicVec()

GSTLEARN_EXPORT VectorDouble ut_sphericalHarmonicVec	(	int	n,
		int	k,
		VectorDouble	theta,
		VectorDouble	phi
	)

ut_vandercorput()

GSTLEARN_EXPORT void ut_vandercorput	(	int	n,
		int	flag_sym,
		int	flag_rot,
		int *	ntri_arg,
		double **	coor_arg
	)

Generate a Van Der Corput list of points in R^3

Parameters

[in]	n	Number of points
[in]	flag_sym	Duplicate the samples by symmetry
[in]	flag_rot	Perform a random rotation
[out]	ntri_arg	Number of points
[out]	coor_arg	Array of point coordinates (Dimension: 3*ntri)

vanDerCorput()

GSTLEARN_EXPORT MatrixRectangular* vanDerCorput	(	int	n,
		int	nd
	)

Function to compute the vector Van der Corput sequence or the Halton sequence

Parameters

n	The number of values to be computed
nd	The dimension of output sequence

Returns

A matrix of dimensions [n,nd] with the sequence values (between 0 and 1)

Note

The dimension nd should be lower or equal to 50.