module igamma;

% Author : Daniel Hobbs , University of Bath, 1995 - 1996
%
%--------------------------------------------------------------------------
%
%  The incomplete gamma function.
%
%  igamma_iter_series(a,x,iter,sum,last_term) - iteratively computes the
%		value of an approximation to an infinite series used in
% 		igamma (for x<=1 or x<a).
%
%  igamma_cont_frac(a,x,iter,iter_max) - iteratively computes the value of
%	the continuous fraction used in igamma (for other values of x).
%
%  igamma_eval(a,x) - returns the value at point x of the
%		incomplete gamma function of order ord.
%
%  The incomplete beta function.
%
%  ibeta_cont_frac(iter,iter_max,a,b,x) - recursively computes
%		the value of the continuous fraction used to
%		approximate to the incomplete beta function.
%
%  ibeta_eval(a,b,x) - returns the value of the incomplete beta
%		function with
%		parameters a and b at point x, by approximating to the
%		incomplete beta function using a continued fraction.
%
%--------------------------------------------------------------------------

% --------------------------- global variables ----------------------------

fluid '(sfiterations);

global '(ibeta_max_iter);

algebraic <<

operator igamma , igamma_eval, ibeta, ibeta_eval;

% Set the maximum number of iterations for the continued fraction
% used in ibeta to be 200.

ibeta_max_iter := 200;

% Set up rule definitions for igamma and ibeta functions.

let
{
 igamma(~a,~x) => igamma_eval(a,x)
        when numberp(a) and numberp(x) and a>0 and x>=0 and lisp !*rounded
};

let
{
 ibeta(~a,~b,~x) => ibeta_eval(a,b,x)
        when numberp(a) and numberp(b) and numberp(x) and lisp !*rounded
             and repart(a)>0 and repart(b)>0 and x>=0 and x<=1
};


% Function igamma_iter_series:     --  cum_gamma_iter        x^i
%				   \			-------------
%				   /			(a+1)...(a+i)
%				   --  i=1
% Uses Battacharjee's method (1970) (computed recursively).

expr procedure igamma_iter_series(a,x,iter,sum,last_term);
begin
 scalar value,this_term;

 if (last_term < 10^-(precision(0)+3)) then
  value := sum
 else
 <<
  this_term := (last_term * x / (a+iter));
  value := igamma_iter_series(a,x,iter+1,sum+this_term,this_term)
 >>;

 return value;
end;


% Function igamma_cont_frac:            1   1-a   1   2-a   2
%				       ---  ---  ---  ---  ---  ...
% 					x +  1 +  x +  1 +  x +
% Recursively computes fraction using
% Abramowitz and Stegun's method (1964).

expr procedure igamma_cont_frac(a,x,iter,iter_max);
begin
 scalar value;

 if (iter>iter_max) then
  value := 0
 else
  value := (iter - a)/
              (1 +      (iter/
                    (x + igamma_cont_frac(a,x,iter + 1,iter_max))));

 return value;
end;


% Function igamma_eval: returns the value at point x of the
% incomplete gamma function with order ord.

expr procedure igamma_eval(a,x);
begin
 scalar arg,frac,last_frac,acc,value;

 % Decide whether to use a series expansion or a continued fraction.
 if (x<=1 or x<a+2) then
 <<
  value := (exp(-x) * x^a) * (1 + igamma_iter_series(a,x,1,0,1)) /
	     gamma(a + 1)
 >>
 else
 <<
  % Set required accuracy to be 3 decimal places more than
  % current precision.
  acc := 10 ^ -(precision(0)+3);
  % Obtain a starting value.
  frac := igamma_cont_frac(a,x,1,1);
  sfiterations := 1;
  % Repeat loop until successive results of continued fraction converge.
  repeat
  <<
   sfiterations := sfiterations + 1;
   last_frac := frac;
   frac := igamma_cont_frac(a,x,1,sfiterations)
  >>
  until (last_frac - frac) < acc;

  arg := exp(-x) * x^a / gamma(a);
  value := 1 - arg / (x + frac)
 >>;

 return value;
end;


% Function ibeta_cont_frac: calculates  1   c(2)  c(3)
%				       ---  ----  ----  ...
%				       1 +  1  +  1  +
% where
%        c(2i) =  - (a + i - 1) (b - i)   *   x
%		 ---------------------------------
%		 (a + 2i - 2) (a + 2i - 1) (1 - x)
% and
%      c(2i+1) =  i (a + b + i - 1)   *   x
%		 -----------------------------
%		 (a + 2i - 1) (a + 2i) (1 - x)

expr procedure ibeta_cont_frac(iter,iter_max,a,b,x);
begin
 scalar value,c_odd,c_even;

 if not (fixp(iter) and fixp(iter_max) and numberp(x)) then
  rederr("ibeta_cont_frac called illegally");

 if (iter>iter_max) then
  value := 0
 else
 <<
  c_even := -(a+iter-1)*(b-iter)*x / ((a+2*iter-2)*(a+2*iter-1)*(1-x));
  c_odd := iter*(a+b+iter-1)*x / ((a+2*iter-1)*(a+2*iter)*(1-x));
  value := c_even /
	       (1 + (c_odd /
		       (1 + ibeta_cont_frac(iter+1,iter_max,a,b,x))))
 >>;

 return value;
end;


% Function ibeta_eval: returns the value of the incomplete beta%
% function with parameters a and b at point x. Method due to Muller (1931).

expr procedure ibeta_eval(a,b,x);
begin
 scalar last_value,value,arg,sfiterations;

 if (x=0 or x=1) then
  value := x
 else
 <<
  % 
  if (repart(a+b)-2)*x > (repart(a)-1) then
   value := 1 - ibeta(b,a,1-x)
  else
  <<
   arg := gamma(a+b) * x^a * (1-x)^(b-1) / (a * gamma(a) * gamma(b));
   % A starting point of 30 levels of continued fraction.
   sfiterations := 30;
   % Starting value that will force calculation a second time at least.
   value := -1;
   repeat
   <<
    last_value := value;
    value := arg * (1/(1 + ibeta_cont_frac(1,sfiterations,a,b,x)));
    sfiterations := sfiterations + 10
   >>
   until (abs(value - last_value) < 10^-(precision(0)+3))
    or sfiterations > ibeta_max_iter;
  >>
 >>;

 % Error condition should not occur, but in case it does...
 if sfiterations > ibeta_max_iter then
 write
 "*** Warning: max iteration limit exceeded; result may not be accurate";

 return value;
end;

>>;

endmodule;

end;

