% EXAMPLE3_64 Computations for Example 3.64.
%
f1='%10.0f %6.0f %12.6f %12.3e %6.0f\n';
f2='        timing=%6.3f\n';
fprintf('\n')
global eps0 Nmax M mc mp iq idelta irout AB Z ab0
eps0=10^9*eps
eps1=10^3*eps
fprintf('\n')
disp('         x      n       H(x)        can        Nfm')
tic
for x=[1 5 10 20 40]
%
% Gauss-Laguerre
%
  Nmax=80;
  ab=r_laguerre(Nmax);
  n=0; gl=1; gl1=0;
  while abs(gl1-gl)>eps0
    gl=gl1; n=n+1;
    if n>Nmax, error('n exceeds Nmax'), end
    qtmax=0; Nfm=0;
    xw=gauss(n,ab);
    gl1=0;
    for k=1:n
      t=xw(k,1);
      if t<1
        term=1; s1=1; s=0; j=1;
        while s1~=s
          s=s1; j=j+1;
          term=-t*term/j;
          s1=s+term;
        end
        e=1/s1;
      else
        e=t/(1-exp(-t));
      end
      [y,Nf]=fHL(t,x,eps1);
      qterm=xw(k,2)*e*y/t;
      if abs(qterm)>qtmax, qtmax=abs(qterm); end
      if Nf>Nfm, Nfm=Nf; end
      gl1=gl1+xw(k,2)*e*y/t;
    end
    can=qtmax/abs(gl1);
  end
  fprintf(f1,x,n,gl1,can,Nfm)
end
time=toc;
fprintf(f2,time)
fprintf('\n')
%
% Rational Gauss-Laguerre
%
Nmax=200;
ab0=r_laguerre(Nmax);
mc=1; mp=0; iq=1; idelta=2; irout=1; AB=[0 Inf];
tic
for x=[1 5 10 20 40]
  if x==40, irout=0; end
  n=0; rgl=1; rgl1=0;
  while abs(rgl1-rgl)>eps0
    rgl=rgl1; n=n+1;
    if n>Nmax, error('n exceeds Nmax'), end
    M=2*floor((n+1)/2);
    for m=1:2:M-1
      Z(m,1)=-1/((m+1)*pi*i); Z(m,2)=1;
      Z(m+1,1)=1/((m+1)*pi*i); Z(m+1,2)=1;
    end
    [abmod,Ncap,kount]=r_mod(n,ab0);
    xw=gauss_rational(n,abmod);
    rgl1=0;
    qtm=0; Nfm=0;
    for k=1:n
      t=xw(k,1);
      if t<1
        term=1; s1=1; s=0; j=1;
        while s1~=s
          s=s1; j=j+1;
          term=-t*term/j;
          s1=s+term;
        end
        e=1/s1;
      else
        e=t/(1-exp(-t));
      end
      [y,Nf]=fHL(t,x,eps1);
      qterm=xw(k,2)*e*y/t;
      if abs(qterm)>qtm, qtm=abs(qterm); end
      if Nf>Nfm, Nfm=Nf; end
      rgl1=rgl1+xw(k,2)*e*y/t;
    end
    can=qtm/abs(rgl1);
  end
  fprintf(f1,x,n,rgl1,can,Nfm)
end
time=toc;
fprintf(f2,time)
fprintf('\n')
%
% Gauss-Einstein
%
Nmax=40;
load -ascii abeinstein1;
ab=abeinstein1(1:Nmax,:);
tic
for x=[1 5 10 20 40]
  n=0; ge=1; ge1=0;
  while abs(ge1-ge)>eps0
    ge=ge1; n=n+1;
    if n>Nmax, error('n exceeds Nmax'), end
    qtmax=0; Nfm=0;
    xw=gauss(n,ab);
    ge1=0;
    for k=1:n
      t=xw(k,1);
      [y,Nf]=fHL(t,x,eps1);
      qterm=xw(k,2)*y/t;
      if abs(qterm)>qtmax, qtmax=abs(qterm); end
      if Nf>Nfm, Nfm=Nf; end
      ge1=ge1+xw(k,2)*y/t;
    end
    can=qtmax/abs(ge1);
  end
  fprintf(f1,x,n,ge1,can,Nfm)
end
time=toc;
fprintf(f2,time)
