% EXAMPLE3_65 Computations for Example 3.65.
%
f1='%10.0f %20.15f %12.4e %12.4e\n';
f2='        timing=%6.3f\n';
fprintf('\n')
Sexact=1-log(2);
%
% Gauss-Laguerre
%
N=50;
disp('         n            S              err          can ')
ab=r_laguerre(N);
tic
for n=1:N
  qtmax=0;
  xw=gauss(n,ab);
  gl=0;
  for k=1:n
    t=xw(k,1);
    qterm=xw(k,2)/(exp(t)+1);
    if abs(qterm)>qtmax, qtmax=abs(qterm); end
    gl=gl+xw(k,2)/(exp(t)+1);
  end
  can=qtmax/abs(gl);
  err=abs((gl-Sexact)/Sexact);
  fprintf(f1,n,gl,err,can)
end
time=toc;
fprintf(f2,time)
fprintf('\n')
%
% Rational Gauss-Laguerre
%
global eps0 Nmax M mc mp iq idelta irout AB Z ab0
eps0=100*eps; Nmax=200;
ab0=r_laguerre(Nmax);
mc=1; mp=0; iq=1; idelta=2; irout=1; AB=[0 Inf];
N=10;
tic
for n=1:N
  M=2*floor((n+1)/2);
  for m=1:2:M-1
    Z(m,1)=-2/((m+1)*pi*i); Z(m,2)=1;
    Z(m+1,1)=2/((m+1)*pi*i); Z(m+1,2)=1;
  end
  [abmod,Ncap,kount]=r_mod(n,ab0);
  xw=gauss_rational(n,abmod);
  rgl=0;
  qtm=0;
  for k=1:n
    t=xw(k,1);
    qterm=xw(k,2)/(exp(t)+1);
    if abs(qterm)>qtm, qtm=abs(qterm); end
    rgl=rgl+xw(k,2)/(exp(t)+1);
  end
  can=qtm/abs(rgl);
  err=abs((rgl-Sexact)/Sexact);
  fprintf(f1,n,rgl,err,can)
end
time=toc;
fprintf(f2,time)
fprintf('\n')
%
% Gauss-Fermi
%
N=20;
load -ascii abfermi1;
ab=abfermi1(1:N,:);
tic
for n=1:N
  qtmax=0;
  xw=gauss(n,ab);
  ge=0;
  for k=1:n
    t=xw(k,1);
    qterm=xw(k,2)*exp(-t);
    if abs(qterm)>qtmax, qtmax=abs(qterm); end
    ge=ge+xw(k,2)*exp(-t);
  end
  can=qtmax/abs(ge);
  err=abs((ge-Sexact)/Sexact);
  fprintf(f1,n,ge,err,can)
end
time=toc;
fprintf(f2,time)