[dpoch1result,a,x]=dpoch1(a,x);

function [dpoch1result,a,x]=dpoch1(a,x);
dpoch1result=[];
persistent absa absx alneps alnvar b bern binv bp first firstCall gbern gbk i ii incr j k ndx nterms pi poly1 q rho sinpx2 sinpxx sqtbig term trig var var2 ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(incr), incr=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(ndx), ndx=0; end;
if isempty(nterms), nterms=0; end;
%***BEGIN PROLOGUE  DPOCH1
%***PURPOSE  Calculate a generalization of Pochhammer's symbol starting
%            from first order.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1, C7A
%***TYPE      doubleprecision (POCH1-S, DPOCH1-D)
%***KEYWORDS  FIRST ORDER, FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate a doubleprecision generalization of Pochhammer's symbol
% for doubleprecision A and X for special situations that require
% especially accurate values when X is small in
%        POCH1(A,X) = (POCH(A,X)-1)/X
%                   = (GAMMA(A+X)/GAMMA(A) - 1.0)/X .
% This specification is particularly suited for stably computing
% expressions such as
%        (GAMMA(A+X)/GAMMA(A) - GAMMA(B+X)/GAMMA(B))/X
%             = POCH1(A,X) - POCH1(B,X)
% Note that POCH1(A,0.0) = PSI(A)
%
% When ABS(X) is so small that substantial cancellation will occur if
% the straightforward formula is used, we use an expansion due
% to Fields and discussed by Y. L. Luke, The Special Functions and Their
% Approximations, Vol. 1, Academic Press, 1969, page 34.
%
% The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
%        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
% In order to maintain significance in POCH1, we write for positive a
%        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
%                       = 1.0 + Q*EXPREL(Q) .
% Likewise the polynomial is written
%        POLY = 1.0 + X*POLY1(A,X) .
% Thus,
%        POCH1(A,X) = (POCH(A,X) - 1) / X
%                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, DCOT, DEXPRL, DPOCH, DPSI, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770801  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)
%***end PROLOGUE  DPOCH1
if isempty(absa), absa=0; end;
if isempty(absx), absx=0; end;
if isempty(alneps), alneps=0; end;
if isempty(alnvar), alnvar=0; end;
if isempty(b), b=0; end;
if isempty(bern), bern=zeros(1,20); end;
if isempty(binv), binv=0; end;
if isempty(bp), bp=0; end;
if isempty(gbern), gbern=zeros(1,21); end;
if isempty(gbk), gbk=0; end;
if isempty(pi), pi=0; end;
if isempty(poly1), poly1=0; end;
if isempty(q), q=0; end;
if isempty(rho), rho=0; end;
if isempty(sinpxx), sinpxx=0; end;
if isempty(sinpx2), sinpx2=0; end;
if isempty(sqtbig), sqtbig=0; end;
if isempty(term), term=0; end;
if isempty(trig), trig=0; end;
if isempty(var), var=0; end;
if isempty(var2), var2=0; end;
if isempty(first), first=false; end;
if firstCall,   bern(1)=[+.833333333333333333333333333333333d-1];  end;
if firstCall,   bern(2)=[-.138888888888888888888888888888888d-2];  end;
if firstCall,   bern(3)=[+.330687830687830687830687830687830d-4];  end;
if firstCall,   bern(4)=[-.826719576719576719576719576719576d-6];  end;
if firstCall,   bern(5)=[+.208767569878680989792100903212014d-7];  end;
if firstCall,   bern(6)=[-.528419013868749318484768220217955d-9];  end;
if firstCall,   bern(7)=[+.133825365306846788328269809751291d-10];  end;
if firstCall,   bern(8)=[-.338968029632258286683019539124944d-12];  end;
if firstCall,   bern(9)=[+.858606205627784456413590545042562d-14];  end;
if firstCall,   bern(10)=[-.217486869855806187304151642386591d-15];  end;
if firstCall,   bern(11)=[+.550900282836022951520265260890225d-17];  end;
if firstCall,   bern(12)=[-.139544646858125233407076862640635d-18];  end;
if firstCall,   bern(13)=[+.353470703962946747169322997780379d-20];  end;
if firstCall,   bern(14)=[-.895351742703754685040261131811274d-22];  end;
if firstCall,   bern(15)=[+.226795245233768306031095073886816d-23];  end;
if firstCall,   bern(16)=[-.574472439520264523834847971943400d-24];  end;
if firstCall,   bern(17)=[+.145517247561486490186626486727132d-26];  end;
if firstCall,   bern(18)=[-.368599494066531017818178247990866d-28];  end;
if firstCall,   bern(19)=[+.933673425709504467203255515278562d-30];  end;
if firstCall,   bern(20)=[-.236502241570062993455963519636983d-31];  end;
if firstCall,   pi=[3.141592653589793238462643383279503d0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DPOCH1
if( first )
sqtbig = 1.0d0./sqrt(24.0d0.*d1mach(1));
alneps = log(d1mach(3));
end;
first = false;
%
if( x==0.0d0 )
[ dpoch1result ,a]=dpsi(a);
end;
if( x==0.0d0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
absx = abs(x);
absa = abs(a);
if( absx>0.1d0.*absa )
%
dpoch1result =(dpoch(a,x)-1.0d0)./x;
elseif( absx.*log(max(absa,2.0d0))>0.1d0 ) ;
dpoch1result =(dpoch(a,x)-1.0d0)./x;
else;
%
bp = a;
if( a<(-0.5d0) )
bp = 1.0d0 - a - x;
end;
incr = 0;
if( bp<10.0d0 )
incr = fix(11.0d0 - bp);
end;
b = bp + incr;
%
var = b + 0.5d0.*(x-1.0d0);
alnvar = log(var);
q = x.*alnvar;
%
poly1 = 0.0d0;
if( var<sqtbig )
var2 =(1.0d0./var).^2;
%
rho = 0.5d0.*(x+1.0d0);
gbern(1) = 1.0d0;
gbern(2) = -rho./12.0d0;
term = var2;
poly1 = gbern(2).*term;
%
nterms = fix(-0.5d0.*alneps./alnvar + 1.0d0);
if( nterms>20 )
xermsg('SLATEC','DPOCH1','NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD',1,2);
end;
if( nterms>=2 )
%
for k = 2 : nterms;
gbk = 0.0d0;
for j = 1 : k;
ndx = fix(k - j + 1);
gbk = gbk + bern(ndx).*gbern(j);
end; j = fix(k+1);
gbern(k+1) = -rho.*gbk./k;
%
term = term.*(2.*k-2-x).*(2.*k-1-x).*var2;
poly1 = poly1 + gbern(k+1).*term;
end; k = fix(nterms+1);
end;
end;
%
poly1 =(x-1.0d0).*poly1;
dpoch1result = dexprl(q).*(alnvar+q.*poly1) + poly1;
%
if( incr~=0 )
%
% WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
% TO OBTAIN DPOCH1(BP,X).
%
for ii = 1 : incr;
i = fix(incr - ii);
binv = 1.0d0./(bp+i);
dpoch1result =(dpoch1result-binv)./(1.0d0+x.*binv);
end; ii = fix(incr+1);
end;
%
if( bp==a )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
% WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A REFLECTION
% FORMULA TO OBTAIN DPOCH1(A,X).
%
sinpxx = sin(pi.*x)./x;
sinpx2 = sin(0.5d0.*pi.*x);
trig = sinpxx.*dcot(pi.*b) - 2.0d0.*sinpx2.*(sinpx2./x);
%
dpoch1result = trig +(1.0d0+x.*trig).*dpoch1result;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DPOCH