/*********************************************************************
Author:  Paul Fackler  paul_fackler@ncsu.edu
         919-515-4535  FACKLER@econbus1.econ.ncsu.edu
Date: Fri, 25 Aug 95
Provided without guarantee for public non-commercial use.
*********************************************************************

           THE LOGGAMMA FUNCTION AND ITS DERIVATIVES

**********************************************************************

  loggamma     (x) log of the gamma function (10 place accuracy)
  psi          (x) psi function (derivative of the log gamma function)
  trigamma     (x) second derivative of the loggamma function
  tetragamma   (x) third derivative of the loggamma function
  pentagamma   (x) fourth derivative of the loggamma function
  polygamma    (x,n) nth derivative of the loggamma function
  psiinv       (p) inverse psi function



Returns the log of the gamma function.  Source: Pike, M.C., and
I.D. Hill, Algorithm 291, Communications of the ACM, 9,9:p.684 (Sept, 1966).
Accepts X>0, accuracy to 10 decimal places.  Uses Sterling's
formula.  Requires workspace at least 5 times the size of X.
*********************************************************************/

proc loggamma(x);
   local z;
     x=x+6;
     z=1/(x.*x);
     z=(((-0.000595238095238*z+0.000793650793651)
          .*z-0.002777777777778).*z+0.083333333333333)./x;
     z=(x-0.5).*ln(x)-x+0.918938533204673+z
           -ln(x-1)-ln(x-2)-ln(x-3)-ln(x-4)-ln(x-5)-ln(x-6);
   retp(z);
endp;

@ Calculates the value of the psi (digamma) function, the first
  derivative of the loggamma function.  Accepts positive matrices.
  From Abromowitz and Stegun.
  Uses formula 6.3.18 with recurrence formula 6.3.5.
  Requires workspace at least 5 times the size of X.
@
proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

@
  3 procs to calculate the value of the trigamma, tetragamma, and
pentagamma
  functions, the second, third, and fourth derivatives of the loggamma
  function. Accepts positive matrices.
  From Abromowitz and Stegun.
  Uses formulas 6.4.11 and 6.4.12 with recurrence formula 6.4.6.
  Each requires workspace at least 5 times the size of X.
@
proc trigamma(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((0.075757575757576*p-0.033333333333333).*p+0.0238095238095238)
     .*p-0.033333333333333).*p+0.166666666666667).*p+1)./x+0.5*p;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=1/(x.*x)+p;
  endo;
  retp(p);
endp;

proc tetragam(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((-0.833333333333333*p+0.3).*p-0.1666666666666667).*p
     +0.1666666666666667).*p-0.5).*p-1).*p;
  p=-1/(x^3)+p;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=-2/(x^3)+p;
  endo;
  retp(p);
endp;

proc pentagam(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((10*p-3).*p+1.3333333333333333).*p-1).*p+2).*p+2+3/x).*p./x;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=6/(x^4)+p;
  endo;
  retp(p);
endp;

@
  Calculates the value of the polygamma function, the Nth derivative
of the
  loggamma function. Accepts positive matrices for X and integer
matrices for
  N. X and N must be ExE conformable.
  From Abromowitz and Stegun.
  Uses formulas 6.4.11 and 6.4.12 with recurrence formula 6.4.6.
@
proc polygam(x,n);
  local p,q,i,b;
  x=x+6;
  p=1/(x.*x);
  q=0;
  let b=0 0.166666666666667 0 -0.033333333333333 0 0.023809523809524 0
-0.03333333333333 0 0.075757575757576 0 -0.253113553113553 0
1.16666666666667;
  i=14; do until i<1.5;
    q=(q+b[i,1]*((i+n-1)!/i!)).*p;
    i=i-2;
  endo;
  p=(-1)^(n-1);
  q=p.*((n-1)! + 0.5*(n!)./x + q)./x^n;
  p=p.*n!;
  i=0; do until i>5;i=i+1;
    x=x-1;
    q=q+p./(x^(n+1));
  endo;
  retp(q);
endp;

@
  Calculates the inverse of the PSI function using the method of
bisection.
  Uses the fact that the exp(p) < psiinv(p) < exp(p)+1.
@
proc psiinv(p);
  local l,x;
  l=1;
  x=exp(p);
  do until l<1e-8;
    x=x+l*sign(p-psi(x));
    l=l/2;
  endo;
  retp(x);
endp;
