[poch1result,a,x]=poch1(a,x);

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

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,9); end;
if isempty(binv), binv=0; end;
if isempty(bp), bp=0; end;
if isempty(gbern), gbern=zeros(1,10); end;
if isempty(gbk), gbk=0; end;
if isempty(pi), pi=0; end;
if isempty(poch1result), poch1result=0; end;
if isempty(poly1), poly1=0; end;
if isempty(q), q=0; end;
if isempty(rho), rho=0; end;
if isempty(sinpx2), sinpx2=0; end;
if isempty(sinpxx), sinpxx=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(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  POCH1
%***PURPOSE  Calculate a generalization of Pochhammer's symbol starting
%            from first order.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1, C7A
%***TYPE      SINGLE PRECISION (POCH1-S, DPOCH1-D)
%***KEYWORDS  FIRST ORDER, FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate a generalization of Pochhammer's symbol 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  COT, EXPREL, POCH, PSI, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770801  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  POCH1
if isempty(first), first=false; end;
if firstCall,   bern(1)=[.83333333333333333e-01];  end;
if firstCall,   bern(2)=[-.13888888888888889e-02];  end;
if firstCall,   bern(3)=[.33068783068783069e-04];  end;
if firstCall,   bern(4)=[-.82671957671957672e-06];  end;
if firstCall,   bern(5)=[.20876756987868099e-07];  end;
if firstCall,   bern(6)=[-.52841901386874932e-09];  end;
if firstCall,   bern(7)=[.13382536530684679e-10];  end;
if firstCall,   bern(8)=[-.33896802963225829e-12];  end;
if firstCall,   bern(9)=[.85860620562778446e-14];  end;
if firstCall,   pi=[3.14159265358979324e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  POCH1
if( first )
sqtbig = 1.0./sqrt(24.0.*r1mach(1));
alneps = log(r1mach(3));
end;
first = false;
%
if( x==0.0 )
[ poch1result ,a]=psi(a);
end;
if( x==0.0 )
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.1.*absa )
%
poch1result =(poch(a,x)-1.0)./x;
elseif( absx.*log(max(absa,2.0))>0.1 ) ;
poch1result =(poch(a,x)-1.0)./x;
else;
%
bp = a;
if( a<(-0.5) )
bp = 1.0 - a - x;
end;
incr = 0;
if( bp<10.0 )
incr = fix(11.0 - bp);
end;
b = bp + incr;
%
var = b + 0.5.*(x-1.0);
alnvar = log(var);
q = x.*alnvar;
%
poly1 = 0.0;
if( var<sqtbig )
var2 =(1.0./var).^2;
%
rho = 0.5.*(x+1.0);
gbern(1) = 1.0;
gbern(2) = -rho./12.0;
term = var2;
poly1 = gbern(2).*term;
%
nterms = fix(-0.5.*alneps./alnvar + 1.0);
if( nterms>9 )
xermsg('SLATEC','POCH1','NTERMS IS TOO BIG, MAYBE R1MACH(3) IS BAD',1,2);
end;
if( nterms>=2 )
%
for k = 2 : nterms;
gbk = 0.0;
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.0).*poly1;
poch1result = exprel(q).*(alnvar+q.*poly1) + poly1;
%
if( incr~=0 )
%
% WE HAVE POCH1(B,X).  BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
% TO OBTAIN POCH1(BP,X).
%
for ii = 1 : incr;
i = fix(incr - ii);
binv = 1.0./(bp+i);
poch1result =(poch1result-binv)./(1.0+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 POCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A REFLECTION
% FORMULA TO OBTAIN POCH1(A,X).
%
sinpxx = sin(pi.*x)./x;
sinpx2 = sin(0.5.*pi.*x);
trig = sinpxx.*cot(pi.*b) - 2.0.*sinpx2.*(sinpx2./x);
%
poch1result = trig +(1.0+x.*trig).*poch1result;
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 POCH