[dpochresult,a,x]=dpoch(a,x);

function [dpochresult,a,x]=dpoch(a,x);
dpochresult=[];
persistent absa absax alnga alngax ax b firstCall i ia n pi sgnga sgngax ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(ia), ia=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  DPOCH
%***PURPOSE  Evaluate a generalization of Pochhammer's symbol.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1, C7A
%***TYPE      doubleprecision (POCH-S, DPOCH-D)
%***KEYWORDS  FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate a doubleprecision generalization of Pochhammer's symbol
% (A)-sub-X = GAMMA(A+X)/GAMMA(A) for doubleprecision A and X.
% For X a non-negative integer, POCH(A,X) is just Pochhammer's symbol.
% This is a preliminary version that does not handle wrong arguments
% properly and may not properly handle the case when the result is
% computed to less than half of doubleprecision.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D9LGMC, DFAC, DGAMMA, DGAMR, DLGAMS, DLNREL, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770701  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  DPOCH
if isempty(absa), absa=0; end;
if isempty(absax), absax=0; end;
if isempty(alnga), alnga=0; end;
if isempty(alngax), alngax=0; end;
if isempty(ax), ax=0; end;
if isempty(b), b=0; end;
if isempty(pi), pi=0; end;
if isempty(sgnga), sgnga=0; end;
if isempty(sgngax), sgngax=0; end;
if firstCall,   pi=[3.141592653589793238462643383279503d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DPOCH
ax = a + x;
if( ax<=0.0d0 )
if( fix(ax)==ax )
%
if( a>0.0D0 || fix(a)~=a )
xermsg('SLATEC','DPOCH','A+X IS NON-POSITIVE INTEGER BUT A IS NOT',2,2);
end;
%
% WE KNOW HERE THAT BOTH A+X AND A ARE NON-POSITIVE INTEGERS.
%
dpochresult = 1.0d0;
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;
%
n = fix(x);
if( min(a+x,a)<(-20.0d0) )
%
dpochresult =(-1.0d0).^n.*exp((a-0.5d0).*dlnrel(x./(a-1.0d0))+x.*log(-a+1.0d0-x)-x+d9lgmc(-a+1.0d0)-d9lgmc(-a-x+1.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;
else;
%
ia = fix(a);
dpochresult =(-1.0d0).^n.*dfac(-ia)./dfac(-ia-n);
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;
end;
end;
%
% A+X IS NOT ZERO OR A NEGATIVE INTEGER.
%
dpochresult = 0.0d0;
if( a<=0.0d0 && fix(a)==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;
%
n = fix(abs(x));
if( (n)~=x || n>20 )
%
absax = abs(a+x);
absa = abs(a);
if( max(absax,absa)<=20.0d0 )
dpochresult = dgamma(a+x).*dgamr(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;
%
elseif( abs(x)>0.5d0.*absa ) ;
%
[dumvar1,alngax,sgngax]=dlgams(a+x,alngax,sgngax);
[a,alnga,sgnga]=dlgams(a,alnga,sgnga);
dpochresult = sgngax.*sgnga.*exp(alngax-alnga);
else;
%
% ABS(X) IS SMALL AND BOTH ABS(A+X) AND ABS(A) ARE LARGE.  THUS,
% A+X AND A MUST HAVE THE SAME SIGN.  FOR NEGATIVE A, WE USE
% GAMMA(A+X)/GAMMA(A) = GAMMA(-A+1)/GAMMA(-A-X+1) *
% SIN(PI*A)/SIN(PI*(A+X))
%
b = a;
if( b<0.0d0 )
b = -a - x + 1.0d0;
end;
dpochresult = exp((b-0.5d0).*dlnrel(x./b)+x.*log(b+x)-x+d9lgmc(b+x)-d9lgmc(b));
if( a<0.0d0 && dpochresult~=0.0d0 )
dpochresult = dpochresult./(cos(pi.*x)+dcot(pi.*a).*sin(pi.*x));
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
return;
end;
else;
%
% X IS A SMALL NON-POSITIVE INTEGER, PRESUMMABLY A COMMON CASE.
%
dpochresult = 1.0d0;
if( n==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;
for i = 1 : n;
dpochresult = dpochresult.*(a+i-1);
end; i = fix(n+1);
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 DPOCO