[pochresult,a,x]=poch(a,x);

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

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(pochresult), pochresult=0; end;
if isempty(sgnga), sgnga=0; end;
if isempty(sgngax), sgngax=0; end;
if isempty(i), i=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  POCH
%***PURPOSE  Evaluate a generalization of Pochhammer's symbol.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C1, C7A
%***TYPE      SINGLE PRECISION (POCH-S, DPOCH-D)
%***KEYWORDS  FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate a generalization of Pochhammer's symbol
% (A)-sub-X = GAMMA(A+X)/GAMMA(A).  For X a non-negative integer,
% POCH(A,X) is just Pochhammer's symbol.  A and X are single precision.
% This is a preliminary version.  Error handling when POCH(A,X) is
% less than half precision is probably incorrect.  Grossly incorrect
% arguments are not handled properly.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  ALGAMS, ALNREL, FAC, GAMMA, GAMR, R9LGMC, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770701  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  POCH
if firstCall,   pi=[3.141592653589793238e0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  POCH
ax = a + x;
if( ax<=0.0 )
if( fix(ax)==ax )
%
if( a>0.0 || fix(a)~=a )
xermsg('SLATEC','POCH','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.
%
pochresult = 1.0;
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;
%
n = fix(x);
if( min(a+x,a)<(-20.0) )
%
pochresult =(-1.0).^n.*exp((a-0.5).*alnrel(x./(a-1.0))+x.*log(-a+1.0-x)-x+r9lgmc(-a+1.)-r9lgmc(-a-x+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;
else;
%
pochresult =(-1.0).^n.*fac(-fix(a))./fac(-fix(a)-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;
%
% HERE WE KNOW A+X IS NOT ZERO OR A NEGATIVE INTEGER.
%
pochresult = 0.0;
if( a<=0.0 && 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( real(n)~=x || n>20 )
%
absax = abs(a+x);
absa = abs(a);
if( max(absax,absa)<=20.0 )
pochresult = gamma(a+x).*gamr(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.5.*absa ) ;
%
[dumvar1,alngax,sgngax]=algams(a+x,alngax,sgngax);
[a,alnga,sgnga]=algams(a,alnga,sgnga);
pochresult = sgngax.*sgnga.*exp(alngax-alnga);
else;
%
% HERE 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.0 )
b = -a - x + 1.0;
end;
pochresult = exp((b-0.5).*alnrel(x./b)+x.*log(b+x)-x+r9lgmc(b+x)-r9lgmc(b));
if( a<0.0 && pochresult~=0.0 )
pochresult = pochresult./(cos(pi.*x)+cot(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.
%
pochresult = 1.0;
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;
pochresult = pochresult.*(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 POIS3D