% EXAMPLE3_63 Computations for Example 3.63.
%
f1='%10.0f %20.15f %12.4e %12.4e\n';
f2='        timing=%6.3f\n';
fprintf('\n')
Texact=1.8600250792211903072;
%
% Generalized Gauss-Laguerre
%
N=20;
disp('         n            T              err          can ')
ab=r_laguerre(N,-1/2);
tic
for n=1:N
  qtmax=0;
  xw=gauss(n,ab);
  gl=0;
  for k=1:n
    x=xw(k,1);
    if x<1
      term=1; s1=1; s=0; j=1;
      while s1~=s
        s=s1; j=j+1;
        term=-x*term/j;
        s1=s+term;
      end
      e=1/s1;
    else
      e=x/(1-exp(-x));
    end
    qterm=xw(k,2)*e*dawson(sqrt(x))/sqrt(x);
    if abs(qterm)>qtmax, qtmax=abs(qterm); end
    gl=gl+xw(k,2)*e*dawson(sqrt(x))/sqrt(x);
  end
  can=qtmax/abs(gl);
  gl=2*gl/sqrt(pi);
  err=abs((gl-Texact)/Texact);
  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=100;
ab0=r_laguerre(Nmax,-1/2);
mc=1; mp=0; iq=1; idelta=2; irout=1; AB=[0 Inf];
N=8;
tic
for n=1:N
  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);
  rgl=0;
  qtm=0;
  for k=1:n
    x=xw(k,1);
    if x<1
      term=1; s1=1; s=0; j=1;
      while s1~=s
        s=s1; j=j+1;
        term=-x*term/j;
        s1=s+term;
      end
      e=1/s1;
    else
      e=x/(1-exp(-x));
    end
    qterm=xw(k,2)*e*dawson(sqrt(x))/sqrt(x);
    if abs(qterm)>qtm, qtm=abs(qterm); end
    rgl=rgl+xw(k,2)*e*dawson(sqrt(x))/sqrt(x);
  end
  can=qtm/abs(rgl);
  rgl=2*rgl/sqrt(pi);
  err=abs((rgl-Texact)/Texact);
  fprintf(f1,n,rgl,err,can)
end
time=toc;
fprintf(f2,time)
fprintf('\n')
%
% Gauss-Einstein
%
N=20;
load -ascii absqm1einstein1;
ab=absqm1einstein1(1:N,:);
tic
for n=1:N
  qtmax=0;
  xw=gauss(n,ab);
  ge=0;
  for k=1:n
    x=xw(k,1);
    qterm=xw(k,2)*dawson(sqrt(x))/sqrt(x);
    if abs(qterm)>qtmax, qtmax=abs(qterm); end
    ge=ge+xw(k,2)*dawson(sqrt(x))/sqrt(x);
  end
  can=qtmax/abs(ge);
  ge=2*ge/sqrt(pi);
  err=abs((ge-Texact)/Texact);
  fprintf(f1,n,ge,err,can)
end
time=toc;
fprintf(f2,time)
