value=hinstm(i,n,s,t,m)

function value=hinstm(i,n,s,t,m)
% Equs. (8) and (9) in Ref.[2]
value=0;
for k=0:n
    value=value+(2*k+1)*gff(2*m-1,k)/gff(2*m+k,k+1)*Pkmi(k,m,i)*Pkmsi(k,m,s,t);
end

function value=Pkmi(k,m,i)
% Gram polynomials (Equ.(14) in Ref.[2])
if k==-1
    value=0;
elseif k==0
    value=1;
else
    value=2*(2*k-1)/k/(2*m-k)*i*Pkmi(k-1,m,i)-(k-1)*(2*m+k-1)/k/(2*m-k)*Pkmi(k-2,m,i);
end

function value=Pkmsi(k,m,s,i)
% s-th derivative of Gram polynomials (Equ. (15) in Ref.[2])
if s==0
    value=Pkmi(k,m,i);
elseif k==-1 || k==0
    value=0;
else
    value=2*(2*k-1)/k/(2*m-k)*(i*Pkmsi(k-1,m,s,i)+s*Pkmsi(k-1,m,s-1,i))-...
        (k-1)*(2*m+k-1)/k/(2*m-k)*Pkmsi(k-2,m,s,i);       
end   

function value=gff(a,b)
% generalized factorial function (Equ.(10) in Ref.[2])
if b==0
    value=1;
else
    value=prod(a-(0:b-1));
end