% I01INFALTAB The integral I_{0,[1,\inf]}(a,b)
%
f0='%8.0f %20.15f\n';
a=1/6;
%b=0;
b=2;
c=a-b+1;
u=6;
N=80; 
ab=r_jacobi01(N);
for n=56
  xw=gauss(n,ab);
  s=0;
  for k=1:n
    x=1+u*xw(k,1)/a;
    w=(-wofy0(-x*exp(-x))/x)^a;
    s=s+xw(k,2)*w*x^(a-b);
  end
  int1=u*s/a;
  [n int1]'
end
return
if c>0
  int2=gamma(c)*(1-gammainc(u+a,c))/a^c;
else
  ab=r_laguerre(N);
  for n=20
    xw=gauss(n,ab);
    s=0;
    for k=1:n
      s=s+xw(k,2)*(1+(u+xw(k,1))/a)^(a-b);
    end
    int2=exp(-u-a)*s/a;
    [n int2]'
  end
end
[int1 int2]'
int=int1+int2;
fprintf(f0,n,int)