function [psiresult,x]=psi(x);
psiresult=[];
persistent apsics aux dxrel first firstCall i n ntapsi ntpsi pi psi psics xbig y ; if isempty(firstCall),firstCall=1;end;
if isempty(apsics), apsics=zeros(1,16); end;
if isempty(aux), aux=0; end;
if isempty(dxrel), dxrel=0; end;
if isempty(pi), pi=0; end;
if isempty(psiresult), psiresult=0; end;
if isempty(psics), psics=zeros(1,23); end;
if isempty(xbig), xbig=0; end;
if isempty(y), y=0; end;
if isempty(i), i=0; end;
if isempty(n), n=0; end;
if isempty(ntapsi), ntapsi=0; end;
if isempty(ntpsi), ntpsi=0; end;
%***BEGIN PROLOGUE PSI
%***PURPOSE Compute the Psi (or Digamma) function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C7C
%***TYPE SINGLE PRECISION (PSI-S, DPSI-D, CPSI-C)
%***KEYWORDS DIGAMMA FUNCTION, FNLIB, PSI FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% PSI(X) calculates the psi (or digamma) function for real argument X.
% PSI(X) is the logarithmic derivative of the gamma function of X.
%
% Series for PSI on the interval 0. to 1.00000D+00
% with weighted error 2.03E-17
% log weighted error 16.69
% significant figures required 16.39
% decimal places required 17.37
%
% Series for APSI on the interval 0. to 2.50000D-01
% with weighted error 5.54E-17
% log weighted error 16.26
% significant figures required 14.42
% decimal places required 16.86
%
%***REFERENCES (NONE)
%***ROUTINES CALLED COT, CSEVL, INITS, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770401 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900727 Added EXTERNAL statement. (WRB)
% 920618 Removed space from variable names. (RWC, WRB)
%***end PROLOGUE PSI
if isempty(first), first=false; end;
if firstCall, psics(1)=[-.038057080835217922e0]; end;
if firstCall, psics(2)=[.49141539302938713e0]; end;
if firstCall, psics(3)=[-.056815747821244730e0]; end;
if firstCall, psics(4)=[.008357821225914313e0]; end;
if firstCall, psics(5)=[-.001333232857994342e0]; end;
if firstCall, psics(6)=[.000220313287069308e0]; end;
if firstCall, psics(7)=[-.000037040238178456e0]; end;
if firstCall, psics(8)=[.000006283793654854e0]; end;
if firstCall, psics(9)=[-.000001071263908506e0]; end;
if firstCall, psics(10)=[.000000183128394654e0]; end;
if firstCall, psics(11)=[-.000000031353509361e0]; end;
if firstCall, psics(12)=[.000000005372808776e0]; end;
if firstCall, psics(13)=[-.000000000921168141e0]; end;
if firstCall, psics(14)=[.000000000157981265e0]; end;
if firstCall, psics(15)=[-.000000000027098646e0]; end;
if firstCall, psics(16)=[.000000000004648722e0]; end;
if firstCall, psics(17)=[-.000000000000797527e0]; end;
if firstCall, psics(18)=[.000000000000136827e0]; end;
if firstCall, psics(19)=[-.000000000000023475e0]; end;
if firstCall, psics(20)=[.000000000000004027e0]; end;
if firstCall, psics(21)=[-.000000000000000691e0]; end;
if firstCall, psics(22)=[.000000000000000118e0]; end;
if firstCall, psics(23)=[-.000000000000000020e0]; end;
if firstCall, apsics(1)=[-.0204749044678185e0]; end;
if firstCall, apsics(2)=[-.0101801271534859e0]; end;
if firstCall, apsics(3)=[.0000559718725387e0]; end;
if firstCall, apsics(4)=[-.0000012917176570e0]; end;
if firstCall, apsics(5)=[.0000000572858606e0]; end;
if firstCall, apsics(6)=[-.0000000038213539e0]; end;
if firstCall, apsics(7)=[.0000000003397434e0]; end;
if firstCall, apsics(8)=[-.0000000000374838e0]; end;
if firstCall, apsics(9)=[.0000000000048990e0]; end;
if firstCall, apsics(10)=[-.0000000000007344e0]; end;
if firstCall, apsics(11)=[.0000000000001233e0]; end;
if firstCall, apsics(12)=[-.0000000000000228e0]; end;
if firstCall, apsics(13)=[.0000000000000045e0]; end;
if firstCall, apsics(14)=[-.0000000000000009e0]; end;
if firstCall, apsics(15)=[.0000000000000002e0]; end;
if firstCall, apsics(16)=[-.0000000000000000e0]; end;
if firstCall, pi=[3.14159265358979324e0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT PSI
if( first )
[ntpsi ,psics]=inits(psics,23,0.1.*r1mach(3));
[ntapsi ,apsics]=inits(apsics,16,0.1.*r1mach(3));
%
xbig = 1.0./sqrt(r1mach(3));
dxrel = sqrt(r1mach(4));
end;
first = false;
%
y = abs(x);
if( y>=2.0 )
%
% PSI(X) FOR ABS(X) .GE. 2.
%
aux = 0.;
if( y<xbig )
[ aux ,dumvar2,apsics,ntapsi]=csevl(8../y.^2-1.,apsics,ntapsi);
end;
if( x<0. )
psiresult = log(abs(x)) - 0.5./x + aux - pi.*cot(pi.*x);
end;
if( x>0. )
psiresult = log(x) - 0.5./x + aux;
end;
else;
%
% PSI(X) FOR -2. .LT. X .LT. 2.
%
n = fix(x);
if( x<0. )
n = fix(n - 1);
end;
y = x - n;
n = fix(n - 1);
[psiresult ,dumvar2,psics,ntpsi]=csevl(2..*y-1.,psics,ntpsi);
if( n==0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
n = fix(-n);
if( x==0. )
xermsg('SLATEC','PSI','X IS 0',2,2);
end;
if( x<0. && x+n-2==0. )
xermsg('SLATEC','PSI','X IS A NEGATIVE INTEGER',3,2);
end;
if( x<(-0.5) && abs((x-fix(x-0.5))./x)<dxrel )
xermsg('SLATEC','PSI','ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
for i = 1 : n;
psiresult = psiresult - 1.0./(x+i-1);
end; i = fix(n+1);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK PSIFN
