Artifact ade160c0b2c07f892537fe0f67ca2550eb0b4e073808703ad8874476e5f26004:
- Executable file
r37/packages/specfn/sfigamma.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5810) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfigamma.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5810) [annotate] [blame] [check-ins using]
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;