function [dpsiresult,x]=dpsi(x);
dpsiresult=[];
persistent apsics aux dxrel first firstCall i n ntapsi ntpsi pi psics xbig y ; if isempty(firstCall),firstCall=1;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 DPSI
%***PURPOSE Compute the Psi (or Digamma) function.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C7C
%***TYPE doubleprecision (PSI-S, DPSI-D, CPSI-C)
%***KEYWORDS DIGAMMA FUNCTION, FNLIB, PSI FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DPSI calculates the doubleprecision Psi (or Digamma) function for
% doubleprecision argument X. PSI(X) is the logarithmic derivative
% of the Gamma function of X.
%
% Series for PSI on the interval 0. to 1.00000E+00
% with weighted error 5.79E-32
% log weighted error 31.24
% significant figures required 30.93
% decimal places required 32.05
%
%
% Series for APSI on the interval 0. to 1.00000E-02
% with weighted error 7.75E-33
% log weighted error 32.11
% significant figures required 28.88
% decimal places required 32.71
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DCOT, DCSEVL, INITDS, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770601 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (WRB)
% 890911 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 name. (RWC, WRB)
%***end PROLOGUE DPSI
if isempty(psics), psics=zeros(1,42); 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(xbig), xbig=0; end;
if isempty(y), y=0; end;
if isempty(first), first=false; end;
if firstCall, psics(1)=[-.38057080835217921520437677667039d-1]; end;
if firstCall, psics(2)=[+.49141539302938712748204699654277d+0]; end;
if firstCall, psics(3)=[-.56815747821244730242892064734081d-1]; end;
if firstCall, psics(4)=[+.83578212259143131362775650747862d-2]; end;
if firstCall, psics(5)=[-.13332328579943425998079274172393d-2]; end;
if firstCall, psics(6)=[+.22031328706930824892872397979521d-3]; end;
if firstCall, psics(7)=[-.37040238178456883592889086949229d-4]; end;
if firstCall, psics(8)=[+.62837936548549898933651418717690d-5]; end;
if firstCall, psics(9)=[-.10712639085061849855283541747074d-5]; end;
if firstCall, psics(10)=[+.18312839465484165805731589810378d-6]; end;
if firstCall, psics(11)=[-.31353509361808509869005779796885d-7]; end;
if firstCall, psics(12)=[+.53728087762007766260471919143615d-8]; end;
if firstCall, psics(13)=[-.92116814159784275717880632624730d-9]; end;
if firstCall, psics(14)=[+.15798126521481822782252884032823d-9]; end;
if firstCall, psics(15)=[-.27098646132380443065440589409707d-10]; end;
if firstCall, psics(16)=[+.46487228599096834872947319529549d-11]; end;
if firstCall, psics(17)=[-.79752725638303689726504797772737d-12]; end;
if firstCall, psics(18)=[+.13682723857476992249251053892838d-12]; end;
if firstCall, psics(19)=[-.23475156060658972717320677980719d-13]; end;
if firstCall, psics(20)=[+.40276307155603541107907925006281d-14]; end;
if firstCall, psics(21)=[-.69102518531179037846547422974771d-15]; end;
if firstCall, psics(22)=[+.11856047138863349552929139525768d-15]; end;
if firstCall, psics(23)=[-.20341689616261559308154210484223d-16]; end;
if firstCall, psics(24)=[+.34900749686463043850374232932351d-17]; end;
if firstCall, psics(25)=[-.59880146934976711003011081393493d-18]; end;
if firstCall, psics(26)=[+.10273801628080588258398005712213d-18]; end;
if firstCall, psics(27)=[-.17627049424561071368359260105386d-19]; end;
if firstCall, psics(28)=[+.30243228018156920457454035490133d-20]; end;
if firstCall, psics(29)=[-.51889168302092313774286088874666d-21]; end;
if firstCall, psics(30)=[+.89027730345845713905005887487999d-22]; end;
if firstCall, psics(31)=[-.15274742899426728392894971904000d-22]; end;
if firstCall, psics(32)=[+.26207314798962083136358318079999d-23]; end;
if firstCall, psics(33)=[-.44964642738220696772598388053333d-24]; end;
if firstCall, psics(34)=[+.77147129596345107028919364266666d-25]; end;
if firstCall, psics(35)=[-.13236354761887702968102638933333d-25]; end;
if firstCall, psics(36)=[+.22709994362408300091277311999999d-26]; end;
if firstCall, psics(37)=[-.38964190215374115954491391999999d-27]; end;
if firstCall, psics(38)=[+.66851981388855302310679893333333d-28]; end;
if firstCall, psics(39)=[-.11469986654920864872529919999999d-28]; end;
if firstCall, psics(40)=[+.19679385886541405920515413333333d-29]; end;
if firstCall, psics(41)=[-.33764488189750979801907200000000d-30]; end;
if firstCall, psics(42)=[+.57930703193214159246677333333333d-31]; end;
if firstCall, apsics(1)=[-.832710791069290760174456932269d-3]; end;
if firstCall, apsics(2)=[-.416251842192739352821627121990d-3]; end;
if firstCall, apsics(3)=[+.103431560978741291174463193961d-6]; end;
if firstCall, apsics(4)=[-.121468184135904152987299556365d-9]; end;
if firstCall, apsics(5)=[+.311369431998356155521240278178d-12]; end;
if firstCall, apsics(6)=[-.136461337193177041776516100945d-14]; end;
if firstCall, apsics(7)=[+.902051751315416565130837974000d-17]; end;
if firstCall, apsics(8)=[-.831542997421591464829933635466d-19]; end;
if firstCall, apsics(9)=[+.101224257073907254188479482666d-20]; end;
if firstCall, apsics(10)=[-.156270249435622507620478933333d-22]; end;
if firstCall, apsics(11)=[+.296542716808903896133226666666d-24]; end;
if firstCall, apsics(12)=[-.674686886765702163741866666666d-26]; end;
if firstCall, apsics(13)=[+.180345311697189904213333333333d-27]; end;
if firstCall, apsics(14)=[-.556901618245983607466666666666d-29]; end;
if firstCall, apsics(15)=[+.195867922607736251733333333333d-30]; end;
if firstCall, apsics(16)=[-.775195892523335680000000000000d-32]; end;
if firstCall, pi=[3.14159265358979323846264338327950d0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DPSI
if( first )
[ntpsi ,psics]=initds(psics,42,0.1.*real(d1mach(3)));
[ntapsi ,apsics]=initds(apsics,16,0.1.*real(d1mach(3)));
%
xbig = 1.0d0./sqrt(d1mach(3));
dxrel = sqrt(d1mach(4));
end;
first = false;
%
y = abs(x);
%
if( y>10.0d0 )
%
% DPSI(X) FOR ABS(X) .GT. 10.0
%
aux = 0.0d0;
if( y<xbig )
[ aux ,dumvar2,apsics,ntapsi]=dcsevl(2.0d0.*(10.0d0./y).^2-1.0d0,apsics,ntapsi);
end;
%
if( x<0.0d0 )
dpsiresult = log(abs(x)) - 0.5d0./x + aux - pi.*dcot(pi.*x);
end;
if( x>0.0d0 )
dpsiresult = log(x) - 0.5d0./x + aux;
end;
else;
%
% DPSI(X) FOR ABS(X) .LE. 2
%
n = fix(x);
if( x<0.0d0 )
n = fix(n - 1);
end;
y = x - n;
n = fix(n - 1);
[dpsiresult ,dumvar2,psics,ntpsi]=dcsevl(2.0d0.*y-1.0d0,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;
%
if( n>0 )
%
% DPSI(X) FOR X .GE. 2.0 AND X .LE. 10.0
%
for i = 1 : n;
dpsiresult = dpsiresult + 1.0d0./(y+i);
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;
else;
%
n = fix(-n);
if( x==0.0D0 )
xermsg('SLATEC','DPSI','X IS 0',2,2);
end;
if( x<0.0D0 && x+n-2==0.0D0 )
xermsg('SLATEC','DPSI','X IS A NEGATIVE INTEGER',3,2);
end;
if( x<(-0.5D0) && abs((x-fix(x-0.5D0))./x)<dxrel )
xermsg('SLATEC','DPSI','ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER',1,1);
end;
%
for i = 1 : n;
dpsiresult = dpsiresult - 1.0d0./(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;
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 DPSIFN
