[x,m,h,ierr]=dhkseq(x,m,h,ierr);

function [x,m,h,ierr]=dhkseq(x,m,h,ierr);
%***BEGIN PROLOGUE  DHKSEQ
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DBSKIN
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (HKSEQ-S, DHKSEQ-D)
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%   DHKSEQ is an adaptation of subroutine DPSIFN described in the
%   reference below.  DHKSEQ generates the sequence
%   H(K,X) = (-X)**(K+1)*(PSI(K,X) PSI(K,X+0.5))/GAMMA(K+1), for
%            K=0,...,M.
%
%***SEE ALSO  DBSKIN
%***REFERENCES  D. E. Amos, A portable Fortran subroutine for
%                 derivatives of the Psi function, Algorithm 610, ACM
%                 Transactions on Mathematical Software 9, 4 (1983),
%                 pp. 494-502.
%***ROUTINES CALLED  D1MACH, I1MACH
%***REVISION HISTORY  (YYMMDD)
%   820601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
%***end PROLOGUE  DHKSEQ
persistent b firstCall fk fln fn fnp gt hrx i j k mx nx r1m5 rln rxsq s slope t tk trm trmh trmr tst u v wdtol xdmy xh xinc xm xmin yint ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(mx), mx=0; end;
if isempty(nx), nx=0; end;
if isempty(gt), gt=0; end;
if isempty(b), b=zeros(1,22); end;
if isempty(fk), fk=0; end;
if isempty(fln), fln=0; end;
if isempty(fn), fn=0; end;
if isempty(fnp), fnp=0; end;
if isempty(hrx), hrx=0; end;
if isempty(rln), rln=0; end;
if isempty(rxsq), rxsq=0; end;
if isempty(r1m5), r1m5=0; end;
if isempty(s), s=0; end;
if isempty(slope), slope=0; end;
if isempty(t), t=0; end;
if isempty(tk), tk=0; end;
if isempty(trm), trm=zeros(1,22); end;
if isempty(trmh), trmh=zeros(1,25); end;
if isempty(trmr), trmr=zeros(1,25); end;
if isempty(tst), tst=0; end;
if isempty(u), u=zeros(1,25); end;
if isempty(v), v=zeros(1,25); end;
if isempty(wdtol), wdtol=0; end;
if isempty(xdmy), xdmy=0; end;
if isempty(xh), xh=0; end;
if isempty(xinc), xinc=0; end;
if isempty(xm), xm=0; end;
if isempty(xmin), xmin=0; end;
if isempty(yint), yint=0; end;
h_shape=size(h);h=reshape(h,1,[]);
%-----------------------------------------------------------------------
%             SCALED BERNOULLI NUMBERS 2.0*B(2K)*(1-2**(-2K))
%-----------------------------------------------------------------------
if firstCall,   b(1) =[1.00000000000000000d+00];  end;
if firstCall,  b(2) =[-5.00000000000000000d-01];  end;
if firstCall,  b(3) =[2.50000000000000000d-01];  end;
if firstCall,  b(4) =[-6.25000000000000000d-02];  end;
if firstCall,  b(5) =[4.68750000000000000d-02];  end;
if firstCall,  b(6) =[-6.64062500000000000d-02];  end;
if firstCall,  b(7) =[1.51367187500000000d-01];  end;
if firstCall,  b(8) =[-5.06103515625000000d-01];  end;
if firstCall, b(9) =[2.33319091796875000d+00];  end;
if firstCall,  b(10) =[-1.41840972900390625d+01];  end;
if firstCall,  b(11) =[1.09941936492919922d+02];  end;
if firstCall,  b(12) =[-1.05824747562408447d+03];  end;
if firstCall,  b(13) =[1.23842434241771698d+04];  end;
if firstCall,  b(14) =[-1.73160495905935764d+05];  end;
if firstCall,  b(15) =[2.85103429084961116d+06];  end;
if firstCall, b(16) =[-5.45964619322445132d+07];  end;
if firstCall,  b(17) =[1.20316174668075304d+09];  end;
if firstCall,  b(18) =[-3.02326315271452307d+10];  end;
if firstCall,  b(19) =[8.59229286072319606d+11];  end;
if firstCall,  b(20) =[-2.74233104097776039d+13];  end;
if firstCall,  b(21) =[9.76664637943633248d+14];  end;
if firstCall, b(22)=[-3.85931586838450360d+16];  end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT  DHKSEQ
ierr = 0;
wdtol = max(d1mach(4),1.0d-18);
fn = m - 1;
fnp = fn + 1.0d0;
%-----------------------------------------------------------------------
%     COMPUTE XMIN
%-----------------------------------------------------------------------
[r1m5 ]=d1mach(5);
rln = r1m5.*i1mach(14);
rln = min(rln,18.06d0);
fln = max(rln,3.0d0) - 3.0d0;
yint = 3.50d0 + 0.40d0.*fln;
slope = 0.21d0 + fln.*(0.0006038d0.*fln+0.008677d0);
xm = yint + slope.*fn;
mx = fix(fix(xm) + 1);
xmin = mx;
%-----------------------------------------------------------------------
%     GENERATE H(M-1,XDMY)*XDMY**(M) BY THE ASYMPTOTIC EXPANSION
%-----------------------------------------------------------------------
xdmy = x;
xinc = 0.0d0;
if( x<xmin )
nx = fix(fix(x));
xinc = xmin - nx;
xdmy = x + xinc;
end;
rxsq = 1.0d0./(xdmy.*xdmy);
hrx = 0.5d0./xdmy;
tst = 0.5d0.*wdtol;
t = fnp.*hrx;
%-----------------------------------------------------------------------
%     INITIALIZE COEFFICIENT ARRAY
%-----------------------------------------------------------------------
s = t.*b(3);
gt=0;
if( abs(s)>=tst )
tk = 2.0d0;
for k = 4 : 22;
t = t.*((tk+fn+1.0d0)./(tk+1.0d0)).*((tk+fn)./(tk+2.0d0)).*rxsq;
trm(k) = t.*b(k);
if( abs(trm(k))<tst )
gt=1;
break;
end;
s = s + trm(k);
tk = tk + 2.0d0;
end;
if(gt==0)
ierr = 2;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end;
end;
h(m) = s + 0.5d0;
gt=0;
if( m~=1 )
%-----------------------------------------------------------------------
%     GENERATE LOWER DERIVATIVES, I.LT.M-1
%-----------------------------------------------------------------------
for i = 2 : m;
fnp = fn;
fn = fn - 1.0d0;
s = fnp.*hrx.*b(3);
if( abs(s)>=tst )
fk = fnp + 3.0d0;
for k = 4 : 22;
trm(k) = trm(k).*fnp./fk;
if( abs(trm(k))<tst )
gt=1;
break;
end;
s = s + trm(k);
fk = fk + 2.0d0;
end;
if(gt==0)
ierr = 2;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end;
end;
mx = fix(m - i + 1);
h(mx) = s + 0.5d0;
end;
end;
if( xinc==0.0d0 )
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end;
%-----------------------------------------------------------------------
%     RECUR BACKWARD FROM XDMY TO X
%-----------------------------------------------------------------------
xh = x + 0.5d0;
s = 0.0d0;
nx = fix(fix(xinc));
for i = 1 : nx;
trmr(i) = x./(x+nx-i);
u(i) = trmr(i);
trmh(i) = x./(xh+nx-i);
v(i) = trmh(i);
s = s + u(i) - v(i);
end; i = fix(nx+1);
mx = fix(nx + 1);
trmr(mx) = x./xdmy;
u(mx) = trmr(mx);
h(1) = h(1).*trmr(mx) + s;
if( m==1 )
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end;
for j = 2 : m;
s = 0.0d0;
for i = 1 : nx;
trmr(i) = trmr(i).*u(i);
trmh(i) = trmh(i).*v(i);
s = s + trmr(i) - trmh(i);
end; i = fix(nx+1);
trmr(mx) = trmr(mx).*u(mx);
h(j) = h(j).*trmr(mx) + s;
end; j = fix(m+1);
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end %subroutine dhkseq
%DECK DHSTRT