%SGAMMASTAR
%
function y=sgammastar(dig,a,x)
syms y y0 t
digits(dig)
eps0=.5*10^(-dig);
a=vpa(a,dig); x=vpa(x,dig);
if rem(subs(a),1)==0 & subs(a)<0
  y=x^(-subs(a));
  return
end
t=vpa(1/gamma(a+1),dig); y0=vpa(0); y=t; k=0;
while subs(abs((y-y0)/y))>eps0
  k=k+1; y0=y; t=vpa(x*t/(a+k),dig); y=y0+t;
end
y=vpa(exp(-x),dig)*y;
