home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / RMATH.ZIP / SPECFUN.DOC < prev    next >
Encoding:
Text File  |  1990-11-01  |  8.8 KB  |  285 lines
  1. Unit SpecFun/SpecF87.  Special functions for mathematical
  2. calculations. Copyright 1990, by J. W. Rider.
  3.  
  4. Like people, all mathematical functions are "special" in their
  5. own way.  The general idea about this particular collection of
  6. "special" function is combinatorics, or different ways of
  7. counting things.
  8.  
  9. References:
  10.  
  11. [HMF] "The Handbook of Mathematical Functions", edited by M.
  12.         Abramowitz and I. A. Stegun, Government Printing Office,
  13.         1970.
  14.  
  15. [NR] "Numerical Recipes", by W. H. Press, B. P. Flannery, S. A.
  16.         Teukolsky, and W. T. Vetterling, Cambridge University
  17.         Press, 1986.
  18.  
  19. Unit dependencies: uses MATH.TPU or MATH87.TPU for float type
  20. definitions, constants and other functions.
  21.  
  22. An "alternative" algorithm is provided for some of the following
  23. functions.  Usually, this is a "brute force" implementation that
  24. is easier to understand than, but numerically inferior to, the
  25. implementation in the unit.
  26.  
  27.  
  28.                  *** Discrete combinatorics ***
  29.  
  30. These functions are generally learned early in a student's
  31. studies, perhaps in algebra, perhaps in probability.  They arise
  32. in numerous instances of numerical applications.
  33.  
  34.  
  35. FACTORIAL (n:word) :xfloat
  36.  
  37.         Returns "n!" or the number of different ways of ordering
  38.         "n" distinct things. the factorial of n = gamma(n+1) The
  39.         "factorial" function attempts to maintain a permanent
  40.         record of its calculations on the heap.
  41.  
  42.     Alternative algorithm (brute force):
  43.  
  44.     function factorial(n:word):xfloat;
  45.     var i:word; f:xfloat;
  46.     begin
  47.        f:=1; { "0!" and "1!" }
  48.        for i:=n downto 2 do f:=f*i;
  49.        factorial:=f;
  50.     end;
  51.  
  52.  
  53.     Alternative algorithm (classic recursive):
  54.  
  55.     function factorial(n:word):xfloat;
  56.     begin
  57.        if n<=1 then
  58.           factorial:=1
  59.        else
  60.           factorial:=factorial(succ(n))*n;
  61.     end;
  62.  
  63.  
  64. LNFACTORIAL (n:word) :xfloat
  65.  
  66.         Returns the natural logarithm of "n!" = ln(factorial(n)).
  67.         The LnFactorial function attempts to preserve its results
  68.         in an array on the heap.
  69.  
  70.         The LN- prefix versions of the factorial and other
  71.         special functions is provided for two reasons:
  72.  
  73.              1) factorials tend to grow large very quickly.
  74.                 Passing the log of the factorial around permits
  75.                 your programs to handle factorials of much larger
  76.                 numbers than the simpler factorial function
  77.                 would,
  78.  
  79.              2) the current implementation naturally calculates
  80.                 the logarithm. Returning the logarithm avoids the
  81.                 redundant step of taking the exponential or
  82.                 anti-logarithm.
  83.  
  84.  
  85. PERMUTATION (n,k:word) :xfloat
  86.  
  87.         Returns the number of permutations (distinct orderings)
  88.         of "n" things taken "k" at a time.  "Permutation(n,n) =
  89.         Factorial(n)".
  90.  
  91.     Alternative algorithm (brute force):
  92.  
  93.     function permutation(n,k:word):xfloat;
  94.     var i:word; p:xfloat;
  95.     begin
  96.        p:=1; { permutation(n,0) }
  97.        for i:=n downto succ(n-k) do p:=p*i;
  98.        permutation:=p;
  99.     end;
  100.  
  101.  
  102.     Alternative algorithm (using factorials):
  103.  
  104.     function permutation(n,k:word):xfloat;
  105.     begin
  106.        permutation:=factorial(n)/factorial(succ(n-k));
  107.     end;
  108.  
  109.  
  110. COMBINATION (n,k:word) :xfloat
  111.  
  112.         Returns the combination (regardless of order) of "n"
  113.         things taken "k" at a time. (NR calls this "bico" or the
  114.         binomial coefficient function.)
  115.  
  116.  
  117.     Alternative algorithm (using factorials):
  118.  
  119.     function combination(n,k:word):xfloat;
  120.     begin
  121.        combination:=factorial(n)/(factorial(k)*factorial(n-k));
  122.     end;
  123.  
  124.  
  125.  
  126.          *** Continous "combinatorics" ***
  127.  
  128. While the discrete combinatorics take integer-like arguments, the
  129. continuous "combinatorics" extend the functions to incorporate
  130. real-like arguments.  They are are called "combinatorics"
  131. because, for integer arguments, the continuous versions return
  132. the same value as do the discrete ones.
  133.  
  134. GAMMA (x:xfloat) :xfloat
  135.  
  136.     Returns the "gamma" function of "x".  The gamma function is a
  137.     continuous version of the factorial function such that
  138.     Factorial(n)=Gamma(n+1).
  139.  
  140.  
  141. LNGAMMA (x:xfloat) :xfloat
  142.  
  143.     Returns the natural logarithm of "gamma(x)" or
  144.     "ln(gamma(x))".  ( NR calls this "gammln".)  The current
  145.     version is based upon that in NR which based its
  146.     implementation on the results of C. Lanczos.  The SpecFun
  147.     implementation extends the NR algorithm to return "reasonable"
  148.     values for x<1. However, there are negative values of "x"
  149.     that would result in a gamma that is not positive.
  150.     Obviously, taking the logarithm of those results would
  151.     normally cause the algorithm to fail. This implementation
  152.     returns the real value of the logarithm instead.  The "gamma"
  153.     algorithm determines when the returned value should be
  154.     negative.
  155.  
  156.  
  157. BETA (x,y:xfloat) :xfloat
  158.  
  159.     Returns the "beta" function of "x" and "y" (order is not
  160.     important).
  161.  
  162.     Alternative algorithm (using "gamma"):
  163.  
  164.     function beta(x,y:xfloat):xfloat;
  165.     begin
  166.        beta:=gamma(x)*gamma(y)/gamma(x+y);
  167.     end;
  168.  
  169.  
  170. LNBETA (x,y:xfloat) :xfloat
  171.  
  172.     Returns natural logarithm of "beta(x,y)", or ln(beta(x,y)).
  173.  
  174.  
  175.          *** "Incomplete" versions of gamma ***
  176.  
  177. Why "incomplete"?  One of the representations of the gamma
  178. function involves taking the integral of an expression from 0 to
  179. infinity.  If you don't go all the way to infinity, the gamma
  180. function is "incomplete".
  181.  
  182. The incomplete gamma function is important because of its
  183. relationship to the gaussian and other probability distributions.
  184.  
  185.  
  186. IGAMMA (a,x:xfloat) :xfloat
  187.  
  188.     Returns *the* "incomplete gamma" function (unfortunately, the
  189.     related functions are confusingly called "incomplete gamma"
  190.     functions also).   IGamma ranges from 0 to gamma(a) as "x"
  191.     ranges from 0 to positive infinity.
  192.  
  193.  
  194. IGAMMAC (a,x:xfloat) :xfloat
  195.  
  196.     Returns the complement of the "incomplete gamma" function, or
  197.     gamma(a)-igamma(a,x).
  198.  
  199.  
  200. IGAMMAP (a,x:xfloat) :xfloat
  201.  
  202.     Returns the degree of completeness of the "incomplete gamma"
  203.     function, or igamma(a,x)/gamma(a).  NR calls this "gammp".
  204.  
  205.  
  206. IGAMMAQ (a,x:xfloat) :xfloat
  207.  
  208.     Returns the complement of IGammaP, or 1-igammap(a,x). NR
  209.     calls this "gammq".
  210.  
  211.  
  212.          *** "Incomplete" versions of beta ***
  213.  
  214. Why "incomplete"?  One of the representations of the beta
  215. function involves taking the integral of an expression from 0 to
  216. 1.  If you don't go all the way to one, the beta function is
  217. "incomplete".
  218.  
  219. Unlike the complete version, the order of the arguments in the
  220. incomplete beta function is important.
  221.  
  222. The incomplete beta function is important because of its
  223. relationship to the binomial and other probability distributions.
  224.  
  225.  
  226. IBETA (a,b,x:xfloat) :xfloat
  227.  
  228.     Returns *the* "incomplete beta" function (related functions
  229.     are also confusingly referred to as "incomplete beta"
  230.     functions). It ranges from 0 to beta(a,b) as x ranges from 0
  231.     to 1.  In HMF, "ibeta(a,b,x)" would be represented as
  232.  
  233.                              B (a,b)
  234.                               x
  235.  
  236.  
  237. There is no "SpecFun.IBetaC" which would be the complement of the
  238. "IBeta" function.  It would be typically defined as
  239. "beta(a,b)-ibeta(a,b,x)" or "ibeta(b,a,1-x)".  However, it does
  240. not normally arise in practice, and there was no separate
  241. algorithm which made it easier to compute under certain
  242. conditions than what was already available for "IBeta".
  243.  
  244.  
  245. IBETAP (a,b,x:xfloat) :xfloat
  246.  
  247.     Returns the degree of completeness of the "incomplete beta"
  248.     function, "ibeta(a,b,x)/beta(a,b)".  NR calls this "betai".
  249.     In HMF, "ibetap(a,b,x)" would be represented as
  250.  
  251.                             I (a,b)
  252.                              x
  253.  
  254.  
  255. IBETAQ (a,b,x:xfloat) :xfloat
  256.  
  257.      Returns the complement of IBetaP, or "1-betap(a,b,x)" which
  258.      is the same as "betap(b,a,1-x)".
  259.  
  260.  
  261.      *** The "error" function and its complement ***
  262.  
  263. While not strictly a combinatoric routine, the "error" function
  264. is closely related to the "gamma" function, so much so that the
  265. SpecFun implementation uses the same internal code for both.
  266. If the argument is limited to between 0 and positive infinity,
  267. the "error" function is a cumulative probability distribution
  268. function of the absolute value of gaussian random deviates.
  269.  
  270.  
  271. ERF (x:xfloat) :xfloat
  272.  
  273.      Returns the "error" function of "x".  This value ranges from
  274.      -1 to 1, as "x" ranges from minus infinity to plus infinity.
  275.  
  276.  
  277. ERFC (x:xfloat) :xfloat
  278.  
  279.      Returns the complement of the "error" function of "x", or
  280.      "1-erf(x)".  ErfC(x) ranges from 2 to 0 as "x" ranges from
  281.      minus infinity to plus infinity. NR provides two routines
  282.      called "erfc" or "erfcc".  The SpecFun implementation is
  283.      closer to the "erfc" implementation.
  284.  
  285.