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.




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

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



[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;


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

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





    Alternative algorithm (classic recursive):


    function factorial(n:word):xfloat;


       if n<=1 then







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



             2) the current implementation naturally calculates

                the logarithm. Returning the logarithm avoids the

                redundant step of taking the exponential or




PERMUTATION (n,k:word) :xfloat


        Returns the number of permutations (distinct orderings)

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



    Alternative algorithm (brute force):


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

    var i:word; p:xfloat;


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

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





    Alternative algorithm (using factorials):


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






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;







         *** 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




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




BETA (x,y:xfloat) :xfloat


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



    Alternative algorithm (using "gamma"):


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






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




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



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.