Unit SpecFun/SpecF87

Unit SpecFun/SpecF87. Special functions for mathematical

Like people, all mathematical functions are "special" in their

own way. The general idea about this particular collection of

"special" function is combinatorics, or different ways of

counting things.

References:

[HMF] "The Handbook of Mathematical Functions", edited by M.

Abramowitz and I. A. Stegun, Government Printing Office,

1970.

[NR] "Numerical Recipes", by W. H. Press, B. P. Flannery, S. A.

Teukolsky, and W. T. Vetterling, Cambridge University

Press, 1986.

Unit dependencies: uses MATH.TPU or MATH87.TPU for float type

definitions, constants and other functions.

An "alternative" algorithm is provided for some of the following

functions. Usually, this is a "brute force" implementation that

is easier to understand than, but numerically inferior to, the

implementation in the unit.

*** Discrete combinatorics ***

These functions are generally learned early in a student's

studies, perhaps in algebra, perhaps in probability. They arise

in numerous instances of numerical applications.

FACTORIAL (n:word) :xfloat

Returns "n!" or the number of different ways of ordering

"n" distinct things. the factorial of n = gamma(n+1) The

"factorial" function attempts to maintain a permanent

record of its calculations on the heap.

Alternative algorithm (brute force):

function factorial(n:word):xfloat;

var i:word; f:xfloat;

begin

f:=1; { "0!" and "1!" }

for i:=n downto 2 do f:=f*i;

factorial:=f;

end;

Alternative algorithm (classic recursive):

function factorial(n:word):xfloat;

begin

if n<=1 then

factorial:=1

else

factorial:=factorial(succ(n))*n;

end;

LNFACTORIAL (n:word) :xfloat

Returns the natural logarithm of "n!" = ln(factorial(n)).

The LnFactorial function attempts to preserve its results

in an array on the heap.

The LN- prefix versions of the factorial and other

special functions is provided for two reasons:

1) factorials tend to grow large very quickly.

Passing the log of the factorial around permits

your programs to handle factorials of much larger

numbers than the simpler factorial function

would,

2) the current implementation naturally calculates

the logarithm. Returning the logarithm avoids the

redundant step of taking the exponential or

anti-logarithm.

PERMUTATION (n,k:word) :xfloat

Returns the number of permutations (distinct orderings)

of "n" things taken "k" at a time. "Permutation(n,n) =

Factorial(n)".

Alternative algorithm (brute force):

function permutation(n,k:word):xfloat;

var i:word; p:xfloat;

begin

p:=1; { permutation(n,0) }

for i:=n downto succ(n-k) do p:=p*i;

permutation:=p;

end;

Alternative algorithm (using factorials):

function permutation(n,k:word):xfloat;

begin

permutation:=factorial(n)/factorial(succ(n-k));

end;

COMBINATION (n,k:word) :xfloat

Returns the combination (regardless of order) of "n"

things taken "k" at a time. (NR calls this "bico" or the

binomial coefficient function.)

Alternative algorithm (using factorials):

function combination(n,k:word):xfloat;

begin

combination:=factorial(n)/(factorial(k)*factorial(n-k));

end;

*** Continous "combinatorics" ***

While the discrete combinatorics take integer-like arguments, the

continuous "combinatorics" extend the functions to incorporate

real-like arguments. They are are called "combinatorics"

because, for integer arguments, the continuous versions return

the same value as do the discrete ones.

GAMMA (x:xfloat) :xfloat

Returns the "gamma" function of "x". The gamma function is a

continuous version of the factorial function such that

Factorial(n)=Gamma(n+1).

LNGAMMA (x:xfloat) :xfloat

Returns the natural logarithm of "gamma(x)" or

"ln(gamma(x))". ( NR calls this "gammln".) The current

version is based upon that in NR which based its

implementation on the results of C. Lanczos. The SpecFun

implementation extends the NR algorithm to return "reasonable"

values for x<1. However, there are negative values of "x"

that would result in a gamma that is not positive.

Obviously, taking the logarithm of those results would

normally cause the algorithm to fail. This implementation

returns the real value of the logarithm instead. The "gamma"

algorithm determines when the returned value should be

negative.

BETA (x,y:xfloat) :xfloat

Returns the "beta" function of "x" and "y" (order is not

important).

Alternative algorithm (using "gamma"):

function beta(x,y:xfloat):xfloat;

begin

beta:=gamma(x)*gamma(y)/gamma(x+y);

end;

LNBETA (x,y:xfloat) :xfloat

Returns natural logarithm of "beta(x,y)", or ln(beta(x,y)).

*** "Incomplete" versions of gamma ***

Why "incomplete"? One of the representations of the gamma

function involves taking the integral of an expression from 0 to

infinity. If you don't go all the way to infinity, the gamma

function is "incomplete".

The incomplete gamma function is important because of its

relationship to the gaussian and other probability distributions.

IGAMMA (a,x:xfloat) :xfloat

Returns *the* "incomplete gamma" function (unfortunately, the

related functions are confusingly called "incomplete gamma"

functions also). IGamma ranges from 0 to gamma(a) as "x"

ranges from 0 to positive infinity.

IGAMMAC (a,x:xfloat) :xfloat

Returns the complement of the "incomplete gamma" function, or

gamma(a)-igamma(a,x).

IGAMMAP (a,x:xfloat) :xfloat

Returns the degree of completeness of the "incomplete gamma"

function, or igamma(a,x)/gamma(a). NR calls this "gammp".

IGAMMAQ (a,x:xfloat) :xfloat

Returns the complement of IGammaP, or 1-igammap(a,x). NR

calls this "gammq".

*** "Incomplete" versions of beta ***

Why "incomplete"? One of the representations of the beta

function involves taking the integral of an expression from 0 to

1. If you don't go all the way to one, the beta function is

"incomplete".

Unlike the complete version, the order of the arguments in the

incomplete beta function is important.

The incomplete beta function is important because of its

relationship to the binomial and other probability distributions.

IBETA (a,b,x:xfloat) :xfloat

Returns *the* "incomplete beta" function (related functions

are also confusingly referred to as "incomplete beta"

functions). It ranges from 0 to beta(a,b) as x ranges from 0

to 1. In HMF, "ibeta(a,b,x)" would be represented as

B (a,b)

There is no "SpecFun.IBetaC" which would be the complement of the

"IBeta" function. It would be typically defined as

"beta(a,b)-ibeta(a,b,x)" or "ibeta(b,a,1-x)". However, it does

not normally arise in practice, and there was no separate

algorithm which made it easier to compute under certain

conditions than what was already available for "IBeta".

IBETAP (a,b,x:xfloat) :xfloat

Returns the degree of completeness of the "incomplete beta"

function, "ibeta(a,b,x)/beta(a,b)". NR calls this "betai".

In HMF, "ibetap(a,b,x)" would be represented as

I (a,b)

IBETAQ (a,b,x:xfloat) :xfloat

Returns the complement of IBetaP, or "1-betap(a,b,x)" which

is the same as "betap(b,a,1-x)".

*** The "error" function and its complement ***

While not strictly a combinatoric routine, the "error" function

is closely related to the "gamma" function, so much so that the

SpecFun implementation uses the same internal code for both.

If the argument is limited to between 0 and positive infinity,

the "error" function is a cumulative probability distribution

function of the absolute value of gaussian random deviates.

ERF (x:xfloat) :xfloat

Returns the "error" function of "x". This value ranges from

-1 to 1, as "x" ranges from minus infinity to plus infinity.

ERFC (x:xfloat) :xfloat

Returns the complement of the "error" function of "x", or

"1-erf(x)". ErfC(x) ranges from 2 to 0 as "x" ranges from

minus infinity to plus infinity. NR provides two routines

called "erfc" or "erfcc". The SpecFun implementation is

closer to the "erfc" implementation.