%
% sI01INFAB Computes the integral I_{0,[1,inf]}(a,b) from its
%    exact integral representation.
%
%    dig = number of decimal digits carried
%     n0 = estimate of the number of quadrature points needed
%      N = maximum number of quadrature points allowed
% 
% The number n0 can be found by running this routine with,
% say, n0=5, and temporarily changing the statement n=n+1
% in the while-loop to n=n+5. Then use for n0 the value n-5,
% where n is the second output variable and restore the 
% statement n=n+1. Select N=50 and dig<32 since the file abv 
% has only 50 rows and a precision of 32 decimal digits.
%
% The routine loads an array abv, assumed to be a copy of
% the appropriate abv-file, and calls on the routine
% sr_jacobi01.m, both of which must be provided.
%
function [I,n]=sI01infab(dig,n0,N,a,b)
if subs(a)<=0
  error('a out of range')
  return
end
ab1=vpaconvert('abv');
ab1=vpa(ab1,dig);
ab2=sr_jacobi01(dig,N);
ab2=vpa(ab2,dig);
e=exp(vpa(1,dig));
err=1; int1=vpa(0); n=n0;
while err>.5*10^(-dig)
  int0=int1;
  n=n+1
  if n>N
    err('n exceeds N')
    return
  end
  xw1=sgauss(dig,n,ab1(1:n,:));
  xw2=sgauss(dig,n,ab2(1:n,:));
  s1=vpa(0); s2=vpa(0);
  for k=1:n
    x1=vpa(xw1(k,1),dig); x2=vpa(xw2(k,1),dig);
    x2=(1-1/e)*x2;
    w1=vpa(xw1(k,2),dig); w2=vpa(xw2(k,2),dig);
    s1=s1+w1*(1-x1)^vpa(-a+b-1);
    s2=s2+w2*(1-x2)^vpa(a-1)* ...
       (-log(1-x2)/x2)^vpa(a-b+1);
  end
  s2=(1-1/e)*s2;
  int1=vpa((-1+a*(s1+s2))/(a-b+1),dig);
% 
% Use the second err-statement if I is expected to
% be substantially larger than 1 in absolute value.
%
%  err=abs(subs(int1-int0));
  err=abs(subs((int1-int0)/int1));
end
I=int1;