digamma & trigamma functions

[Christoph Mecklenbraeuker]

Hi Guys,

A few days ago, I programmed the 'Digamma-' (aka 'Psi-') and 'Trigamma-'
functions which are part of the more general family of 'Polygamma-Functions'.
{which are recursively defined as derivatives of log(gamma(x))}

The code is tested  both for octave-1.1.0 and Matlab-4.2 on SunOS platform.
(sorry: I'm not up-to-date with the latest version of octave ...)

Caveat: Accuracy of the evaluation for *any* real argument was not my main
concern, because I need the functions only for integer and half-integer 
arguments,
where *exact* finite recursions are readily availlable. So, I only
guarantee full precision for numbers of the form 1/2, 2/2, 3/2, 4/2,
5/2, 6/2, 7/2, ... but 3-4 decimals should be correct anywhere on the
real line. Maybe in near future, I will post an update with increased
accuracy. 

So, here you get my alpha-version of both functions: suggestions for 
improvements are most welcome. I tried to make the code as time
efficient as possible, but a lot of work still needs to be done.

Have fun!
Christoph

--(cut-here)--------------------file:digamma.m----------------------
function psi = digamma(z,method,debug)
%
%  psi = digamma(z)   ... Digamma-Function for real argument z.
%
%  digamma(z) = d/dz log(gamma(z)) = gamma'(z)/gamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for 
%  each element. Results may be inaccurate for real arguments < 10 
%  which are neither integers nor half-integers.
%
%  psi = digamma(z,method)
%
%  possible values for optional argument 'method': 
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%
%  see also:  trigamma, gamma, gammaln, gammainc, specfun

%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: address@hidden), July 1, 1995.


dim = size(z);                          % save original matrix dimension
z = reshape(z,dim(1)*dim(2),1);         % make a column vector
I1 = ones(length(z),1);                 % auxiliary vector of ones

if(nargin==1) 
  method=4; debug=0; 
elseif(nargin==2) 
  debug=0; 
end;

if(debug == 1)                          % if debug==1: track recursion
  [m,n] = size(z);
  fprintf(1,'digamma: method = %d, size(z)=[%d %d],\t min(z)=%f, 
max(z)=%f\n',...
      method,m,n,min(min(z)),max(max(z)));
end;


if(method==1)                           % use 8th order asymptotic expansion 
  if(any(z<1))
    fprintf(1,'Warning: some elements in argument of "digamma(z,1)" are < 1\n');
    fprintf(1,'minimal argument = %g: digamma-result is 
inaccurate!\n',min(min(z)));
  end
  % calculate powers of 1/z :
  w1 = 1./z; w2 = w1.*w1; w4 = w2.*w2; w6 = w2.*w4; w8 = w4.*w4;
  % generate coefficients of expansion: matrix with constant columns
  a = [ -I1/2   -I1/12   I1/120  -I1/252  I1/240 ];
  % make vector of powers of 1/z:
  w = [  w1      w2       w4       w6       w8   ];
  % calculate expansion by summing the ROWS of (a .* w) :
  psi = log(z) + sum((a.*w).').';       
elseif(method==2)
  zmax = max(max(floor(z)));
  psitab = zeros(zmax,1);
  psitab(1) = -0.5772156649015328606;
  for n=1:zmax-1;
    psitab(n+1) = psitab(n) + 1/n;      % generate lookup table
  end;
  psi = psitab(z);
elseif(method==3)
  zmax = max(max(floor(z)));
  psitab = zeros(zmax+1,1);
  psitab(1) = -0.5772156649015328606 - 2*log(2);  % = psi(1/2)
  for n=1:zmax;
    psitab(n+1) = psitab(n) + 2/(2*n-1); % generate lookup table
  end;
  psi = psitab(z+0.5);
elseif(method==4)                       % decide here which method to use
  Less0 = find(z<0);                    % negative arguments evaluated by 
reflexion formula
  Less1 = find(z>0 & z<1);              % values between 0 and 1.
  fraction = rem(z,1);                  % fractional part of arguments
  f2 = rem(2*fraction,1);
  Integers = find(fraction==0 & z>0);   % Index set of positive integer 
arguments
  NegInts  = find(fraction==0 & z<=0);  % Index set of positive integer 
arguments
  HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive 
half-integers
  Reals    = find(f2>1e-7 & z>1);       % Index set of all other arguments > 1
  if(~isempty(Reals)) psi(Reals)    = digamma(z(Reals),1,debug);    end;
  if(~isempty(Less1)) psi(Less1)    = digamma(z(Less1)+2,1,debug) - ...
        1./z(Less1)-1./(z(Less1)+1);end;
  % reflexion formula:
  if(~isempty(Less0)) psi(Less0) = digamma(1-z(Less0),1,debug) - 
pi./tan(pi*z(Less0)); end;
  if(~isempty(Integers)) psi(Integers) = digamma(z(Integers),2,debug); end;
  if(~isempty(HalfInts)) psi(HalfInts) = digamma(z(HalfInts),3,debug); end;
  if(~isempty(NegInts))  psi(NegInts)  = Inf * NegInts; end;
end

psi = reshape(psi,dim(1),dim(2));

return;
--(cut-here)-----------------end-of-file:digamma.m----------------------

and here is the second one:
--(cut-here)--------------------file:trigamma.m----------------------
function y = trigamma(z,method,debug)

%  y = trigamma(z)   ... Trigamma-Function for real positive z
%
%  trigamma(z) = (d/dz)^2 log(gamma(z)) = d/dz digamma(z)
%
%  if 'z' is a matrix, then the digamma-function is evaluated for 
%  each element. Results are inaccurate for real arguments < 10 which are 
%  neither integers nor half-integers.
%
%  y = trigamma(z,method)
%
%  possible values for optional argument 'method':
%    method = 1           : quick asymptotic series expansion (approximate)
%    method = 2           : finite recursion for integer values (exact)
%    method = 3           : finite recursion for half-integer values (exact)
%    method = 4 (default) : automatic selection of 1,2 or 3 for individual
%                           elements in z whichever is appropriate.
%
%  see also: digamma, gamma, gammaln, gammainc, specfun


%  reference: Abramowitz & Stegun, "Handbook of Mathematical Functions"
%             Chapter "Gamma Function and Related Functions" :
%  implemented by: Christoph Mecklenbraeuker
%             (email: address@hidden), July 4, 1995.


dim = size(z);                          % save original matrix dimension
z = reshape(z,dim(1)*dim(2),1);         % make a column vector
I1 = ones(length(z),1);                 % auxiliary vector of ones

if(nargin==1)
  method=4; debug=0;
elseif(nargin==2)
  debug=0;
end;


if(debug == 1)                          % if debug==1: track recursion
  [m,n] =size(z);
  fprintf(1,'trigamma: method = %d, size(z)=[%d %d],\t min(z)=%f, 
max(z)=%f\n',...
      method,m,n,min(min(z)),max(max(z)));
end;

if(method==1)                           % use 9th order asymptotic expansion 
  if(any(z<1))
    fprintf(1,'Warning: some elements in argument of "trigamma(z,1)" are < 
1\n');
    fprintf(1,'minimal argument = %g: trigamma-result is 
inaccurate!\n',min(min(z)));
  end

  % calculate powers of 1/z :
  w1 = 1./z; w2 = w1.*w1; w3 = w1.*w2; w5 = w2.*w3; w7 = w2.*w5; w9 = w2.*w7;
  % generate coefficients of expansion: matrix with constant columns
  a = [ I1   I1/2   I1/6  -I1/30  I1/42 -I1/30];
  % make vector of powers of 1/z:
  w = [ w1   w2     w3     w5      w7    w9];
  % calculate expansion by summing the ROWS of (a .* w) :
  y = sum((a.*w).').';  
elseif(method==2)
  zmax = max(max(floor(z)));
  ytab = zeros(zmax,1);
  ytab(1) = pi^2/6;                     % = psi'(1)
  for n=1:zmax-1;
    ytab(n+1) = ytab(n) - 1/n^2;        % generate lookup table
  end;
  y = ytab(z);
elseif(method==3)
  zmax = max(max(floor(z)));
  ytab = zeros(zmax+1,1);
  ytab(1) = pi^2/2;                     % = psi'(1/2)
  for n=1:zmax;
    ytab(n+1) = ytab(n) - 4/(2*n-1)^2; % generate lookup table
  end;
  y = ytab(z+0.5);
elseif(method==4)                       % decide here which method to use
  Less0 = find(z<0);                    % negative arguments evaluated by 
reflexion formula
  Less1 = find(z>0 & z<1);              % values between 0 and 1.
  fraction = rem(z,1);                  % fractional part of arguments
  f2 = rem(2*fraction,1);
  Integers = find(fraction==0 & z>0);   % Index set of positive integer 
arguments
  NegInts  = find(fraction==0 & z<=0);  % Index set of positive integer 
arguments
  HalfInts = find(abs(fraction-0.5)<1e-7 & z>0); % Index set of positive 
half-integers
  Reals    = find(f2>1e-7 & z>1);       % Index set of all other arguments > 1
  if(~isempty(Reals)) y(Reals)    = trigamma(z(Reals),1,debug);    end;
  if(~isempty(Less1)) y(Less1)    = trigamma(z(Less1)+2,1,debug) + ...
        1./z(Less1).^2+1./(z(Less1)+1).^2;end;
  % reflexion formula:
  if(~isempty(Less0)) y(Less0)= 
-trigamma(1-z(Less0),1,debug)+(pi./sin(pi*z(Less0))).^2; end;
  % integers:
  if(~isempty(Integers)) y(Integers) = trigamma(z(Integers),2,debug); end;
  % half-integers:
  if(~isempty(HalfInts)) y(HalfInts) = trigamma(z(HalfInts),3,debug); end;
  % negative integers:
  if(~isempty(NegInts))  y(NegInts)  = Inf * NegInts; end;
end

  y = reshape(y,dim(1),dim(2));
return;
--(cut-here)--------------------end-of-file:trigamma.m----------------------

-- 
Christoph Mecklenbraeuker             | send electronic mail to:
Lehrstuhl f. Signaltheorie (IC 5/35)  | address@hidden
Ruhr Universitaet Bochum              | Tel: +49(234) 700 6119
D-44780 Bochum, Germany               | Fax: +49(234) 709 4261 
---