% CHALLEPS The SIAM 100-digit challenge integral using 
%    the epsilon algorithm
%
f1='%10.0f %23.15f\n';
tic
n=21; M1=20; M2=13;
tol=.5e-15;
s=zeros(n,1); sa=zeros((n+1)/2,1);
ab0=r_jacobi01(M1);
ab=r_jacobi(M2);
U=zeros(M1,1);
xw=gauss(M1,ab0);
xw=pi*xw/2;  
for i=1:M1
    U(i)=uofx(xw(i,1));
end
s(1)=sum(xw(:,2).*cos(xw(:,1))./(xw(:,1)+U));
sgn=1;
for k=1:n
  if k>=2, s(k)=s(k-1)+sgn*a; end
  U=zeros(M2,1);
  xw=gauss(M2,ab);
  xw=pi*xw/2;
  for i=1:M2
    U(i)=uofx(xw(i,1)+k*pi);
  end
  a=sum(xw(:,2).*cos(xw(:,1))./(xw(:,1)+U+k*pi));
  sgn=-sgn;
end
E=epsalg(n,s);
for i=1:(n+1)/2
  sa(i)=E(n-2*i+2,2*i);
end
time=toc;
fprintf('\n')
disp('         k            sa(k)')
for k=1:(n+1)/2
  fprintf(f1,k,sa(k))
end
fprintf('\n')
fprintf('        time=%6.4f',time)