% TABLE2_10  Computations for Table 2.10.
%
f2='%10.0f   <=%3.0f     <=1e-14       <=1e-14\n';
f3='%17.0f %13.3e %13.3e\n';
f4='                  %13.3e %13.3e\n';
f5='                      %6.3f';
fprintf('\n')
disp('   Errors in the recursion coefficients of discrete')
disp('   Fejer polynomials computed by the Stieltjes')
disp('   procedure')
fprintf('\n')
disp('        N      n        erra          errb')
fprintf('\n')
tic
for N=[40 80 160 320]
  xw=fejer(N);
  abex=zeros(N,2);
  if N==40
    load -ascii bfej40;
    abex(:,2)=bfej40;
  elseif N==80
    load -ascii bfej80;
    abex(:,2)=bfej80;
  elseif N==160
    load -ascii bfej160;
    abex(:,2)=bfej160;
  else
    load -ascii bfej320;
    abex(:,2)=bfej320;
  end
  ab=stieltjes(N,xw);
  errab=zeros(N,2);
  errab(:,1)=abs(ab(:,1)-abex(:,1));
  errab(:,2)=abs(ab(:,2)-abex(:,2));
  erram=max(errab(:,1)); errbm=max(errab(:,2));
  f=find(errab(:,1)>1e-14 | errab(:,2)>1e-14);
  if isempty(f)
    fprintf(f2,N,N)
    fprintf(f4,[erram errbm]')
  else
    errab=errab(f,:);
    fprintf(f2,N,f(1)-2)
    if N==40
      fprintf(f3,[f-1 errab]')
      fprintf(f4,[erram errbm]')
    elseif N==80
      fprintf(f3,[f(1:4:13)-1 errab(1:4:13,:)]')
      fprintf(f4,[erram errbm]')
    elseif N==160
      fprintf(f3,[f(1:7:22)-1 errab(1:7:22,:)]')
      fprintf(f4,[erram errbm]')
    else
      fprintf(f3,[f(1:2:5)-1 errab(1:2:5,:)]')
      fprintf(f4,[erram errbm]')
    end
  end
end
t=toc;
fprintf('\n')
disp('                      timing')
fprintf(f5,t)
