function [cpsiresult,zin]=cpsi(zin);
cpsiresult=[];
persistent bern bound cabsz corr dxrel first firstCall i n ndx nterm pi rbig rmin x y z z2inv ; if isempty(firstCall),firstCall=1;end;
;
if isempty(bern), bern=zeros(1,13); end;
if isempty(bound), bound=0; end;
if isempty(cabsz), cabsz=0; end;
if isempty(dxrel), dxrel=0; end;
if isempty(pi), pi=0; end;
if isempty(rbig), rbig=0; end;
if isempty(rmin), rmin=0; end;
if isempty(x), x=0; end;
if isempty(y), y=0; end;
if isempty(i), i=0; end;
if isempty(n), n=0; end;
if isempty(ndx), ndx=0; end;
if isempty(nterm), nterm=0; end;
%***BEGIN PROLOGUE CPSI
%***PURPOSE Compute the Psi (or Digamma) function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C7C
%***TYPE COMPLEX (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 of X. PSI(X)
% is the logarithmic derivative of the gamma function of X.
%
%***REFERENCES (NONE)
%***ROUTINES CALLED CCOT, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 780501 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)
%***end PROLOGUE CPSI
if isempty(z), z=0; end;
if isempty(z2inv), z2inv=0; end;
if isempty(corr), corr=0; end;
if isempty(first), first=false; end;
if firstCall, bern(1)=[.83333333333333333e-1]; end;
if firstCall, bern(2)=[-.83333333333333333e-2]; end;
if firstCall, bern(3)=[.39682539682539683e-2]; end;
if firstCall, bern(4)=[-.41666666666666667e-2]; end;
if firstCall, bern(5)=[.75757575757575758e-2]; end;
if firstCall, bern(6)=[-.21092796092796093e-1]; end;
if firstCall, bern(7)=[.83333333333333333e-1]; end;
if firstCall, bern(8)=[-.44325980392156863e0]; end;
if firstCall, bern(9)=[.30539543302701197e1]; end;
if firstCall, bern(10)=[-.26456212121212121e2]; end;
if firstCall, bern(11)=[.28146014492753623e3]; end;
if firstCall, bern(12)=[-.34548853937728938e4]; end;
if firstCall, bern(13)=[.54827583333333333e5]; end;
if firstCall, pi=[3.141592653589793e0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT CPSI
if( first )
nterm = fix(-0.30.*log(r1mach(3)));
% MAYBE BOUND = N*(0.1*EPS)**(-1/(2*N-1)) / (PI*EXP(1))
bound = 0.1171.*nterm.*(0.1.*r1mach(3)).^(-1.0./(2.*nterm-1));
dxrel = sqrt(r1mach(4));
rmin = exp(max(log(r1mach(1)),-log(r1mach(2)))+0.011);
rbig = 1.0./r1mach(3);
end;
first = false;
%
z = zin;
x = real(z);
y = imag(z);
if( y<0.0 )
z = conj(z);
end;
%
corr =complex(0.0,0.0);
cabsz = abs(z);
if( x<0.0 || cabsz<=bound )
if( x>=0.0 || abs(y)<=bound )
%
if( cabsz<bound )
%
% use THE RECURSION RELATION FOR ABS(Z) SMALL.
%
if( cabsz<rmin )
xermsg('SLATEC','CPSI','CPSI CALLED WITH Z SO NEAR 0 THAT CPSI OVERFLOWS',2,2);
end;
%
if( x<(-0.5) && abs(y)<=dxrel )
if( abs((z-fix(x-0.5))./x)<dxrel )
xermsg('SLATEC','CPSI','ANSWER LT HALF PRECISION BECAUSE Z TOO NEAR NEGATIVE INTEGER',1,1);
end;
if( y==0.0 && x==fix(x) )
xermsg('SLATEC','CPSI','Z IS A NEGATIVE INTEGER',3,2);
end;
end;
%
n = fix(sqrt(bound.^2-y.^2) - x + 1.0);
for i = 1 : n;
corr = corr - 1.0./z;
z = z + 1.0;
end; i = fix(n+1);
else;
%
% use THE REFLECTION FORMULA FOR REAL(Z) NEGATIVE, ABS(Z) LARGE, AND
% ABS(AIMAG(Y)) SMALL.
%
corr = -pi.*ccot(pi.*z);
z = 1.0 - z;
end;
end;
end;
%
% NOW EVALUATE THE ASYMPTOTIC SERIES FOR SUITABLY LARGE Z.
%
if( cabsz>rbig )
cpsiresult = log(z) + corr;
end;
if( cabsz<=rbig )
%
cpsiresult =complex(0.0,0.0);
z2inv = 1.0./z.^2;
for i = 1 : nterm;
ndx = fix(nterm + 1 - i);
cpsiresult = bern(ndx) + z2inv.*cpsiresult;
end; i = fix(nterm+1);
cpsiresult = log(z) - 0.5./z - cpsiresult.*z2inv + corr;
end;
%
if( y<0.0 )
cpsiresult = conj(cpsiresult);
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',zin); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK CPTSL
