%STURAN_LAGUERRE_EXPLORE Symbolic Gauss-Turan quadrasture formula
%  for Laguerre weight function. Preliminary routine for purposes
%  of exploration.
%
global  s
format long e
hom=0;
s=4; N=16; 
dig=72; eps0=.5e-70; eps1=.5e-14;
fprintf('digits=%3.0f,  s=%2.0f,  N=%3.0f\n',dig,s,N)
%
% select values of n for which complete reults 
% are to be printed
nprint=0;
%
R=vpa(zeros(N,2*N-1)); al0=vpa(zeros(5,1));
x0=[1.596071637983321;2.180607124035125;2.759002709962271; ...
3.333551485269048;3.905451757616910];
for j=1:5
  al0(j)=sroot(dig,@p,@pd,eps0,x0); 
end
R(1,1)=al0(s);
rho(1)=R(1,1);
sab0=sr_laguerre(dig,(s+1)*N);
fprintf('eps1=%12.4e\n',eps1)
for n=1:N
  print=nnz(nprint-n)<size(nprint,1);
  if n==1
    rho_init=rho(1);
  elseif n==2
    rho_init=[rho(1);vpa(linfunGTL(2));vpa(quadfunGTL(2))];
%
%   for the initial run (with N=4)
%
%    rho_init=[rho(1);3.7*sab0(2,1);11.3*sab0(2,2)];
%
  else
    rho_init(2*n-1)=vpa(quadfunGTL(n));
%
%   for the initial run (with N=4)
%
%    if n==3, rho_init(2*n-1)=14.4*sab0(3,2); end
%    if n==4, rho_init(2*n-1)=15.3*sab0(4,2); end
%
    for k=n+1:2*n-2
      rho_init(k)=rho(k-1);
    end
    rho_init(n)=vpa(linfunGTL(n));
%   for the initial run (with N=4)
%
%    if n==3, rho_init(n)=3.9*sab0(3,1); end
%    if n==4, rho_init(n)=3.9*sab0(4,1); end
%
    for k=1:n-1
      rho_init(k)=rho(k);
    end
  end
  fprintf('\n')
  tic
  [sxw,rho,it,Jmax,condJmax]=sturan(dig,n,s,eps1,sab0,rho_init,hom,print);
  time=toc;
  fprintf('n=%2.0f\n',n)
  fprintf('time=%8.2f\n',time)
  fprintf('it=%2.0f\n',it)
  fprintf('Jmax=%12.4e\n',Jmax)
  fprintf('condJmax=%12.4e\n',condJmax)
  eps0=.5*10^(-(dig-log10(condJmax)-2));
  if eps0>eps1
    error('eps1 too small to cope with ill-conditioning')
  end
  R(n,1:2*n-1)=rho';
%
% the constant terms in the linear resp. quadratic
% function determining initial values in Newton's method
%
%  if n==2 & N==4
%    fprintf('R(2,2)=%8.4f\n',subs(R(2,2)))
%    fprintf('R(2,3)=%8.4f\n',subs(R(2,3)))
%  end
%
  if print
    disp('                       alpha')   
    subs(R(n,1:n))'
    fprintf('\n')
    disp('                        beta')
    subs([0;R(n,n+1:2*n-1)'])
    fprintf('\n')
    for sig=1:2*s+2
      if sig==1
        disp('                       nodes')
      else
        disp('                      weights')
      end
      k=1:n;
      subs(sxw(:,sig))
      fprintf('\n')
    end
  end
end
%
% when N=4 determine linear and quadratic functions for
% computing initial approximations in Newton's method
%
%if N==4
%  a1=(R(4,4)-R(2,2))/2
%  b1=R(3,5)-R(2,3), d=R(4,7)-R(3,5); b2=(d-b1)/2
%
% the linear resp. quadratic functions are
%     R(2,2)+a1*(n-2),
%     R(2,3)+b1*(n-2)+b2*(n-2)*(n-3)
%
%end
%n=(1:N-1)';
%n=(1:N-6)';
%for k=1:N-1
%  for m=2:N+3-k
%    Rp(m-1,k)=subs(R(m+k-1,2*m+k-2));
%  end
C=zeros(N,5);
merra=zeros(N-1,1);
for k=0:N-1
  if k<N-1
    d=diag(subs(R),-k); sd=size(d,1);
    dd=(d(sd)-d(1))/(sd-1); C(k+1,1)=dd; C(k+1,2)=d(1)-dd;
    m=(1:sd)'; x=linspace(1,N)';
    merra(k+1)=max(abs((C(k+1,1)*m+C(k+1,2)-d)./d));
    if k==0
      figure
      plot(m,d,'*')
      hold on
      plot(x,C(k+1,1)*x+C(k+1,2))
      hold off
    end
  else
    C(k+1,1)=0; C(k+1,2)=subs(R(k+1,1));
  end
end
merrb=zeros(N-3,1);
for k=1:N-1
  d=zeros(N-k,1);
  for m=k+1:N
    d(m-k)=R(m,2*m-k);
  end
  if k<N-2
    mc=ceil((N-k+1)/2);
    d1=(d(mc)-d(1))/(mc-1); d2=(d(N-k)-d(mc))/(N-k-mc);
    dd=(d2-d1)/(N-k-1);
    C(k,3)=dd; C(k,4)=d1-dd*(mc+1); C(k,5)=d(1)+dd*mc-d1;
    m=(1:N-k)';
    merrb(k)=max(abs((C(k,3)*m.^2+C(k,4)*m+C(k,5)-d)./d));
    if k==1
      figure
      plot(m,d,'*')
      hold on
      plot(x,C(k,3)*x.^2+C(k,4)*x+C(k,5))
      hold off
    end
  elseif k==N-2
    C(k,3)=0; C(k,4)=d(2)-d(1); C(k,5)=2*d(1)-d(2);
    m=(1:2)';
%    figure
%    plot(m,d,'*')
%    hold on 
%    plot(x,C(k,4)*x+C(k,5))
%    hold off
  else
    C(k,3)=0; C(k,4)=0; C(k,5)=R(N,N+1);
  end
end
fprintf('n=%2.0f\n',n)
merra
merrb
subs(R)
C
