Unit
SpecFun/SpecF87. Special functions for
mathematical

calculations.
Copyright 1990, by J. W. Rider.

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)

x

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)

x

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.