function r=NewtonCorrection(x,alpha,beta); 
% NEWTONCORRECTION computes directly the Newton correction
%    r=NewtonCorrection(x,alpha,beta) evaluates the ratio of the 
%    polynomial and its derivative given by the three-term 
%    recurrence p_i(x)=(x-alpha_i)p_{i-1}(x)-beta_{i-1}^2
%    p_{i-2}(x) for the Newton correction

n=length(alpha); 
p1s=0; p1=1; p2s=1; p2=x-alpha(1);
for k=1:n-1 
  p0s=p1s; p0=p1; p1s=p2s; p1=p2; 
  p2s=p1+(x-alpha(k+1))*p1s-beta(k)^2*p0s; 
  p2=(x-alpha(k+1))*p1-beta(k)^2*p0;
  maxx=abs(p2)+abs(p2s); 
  if maxx>1e20; d=1e-20; 
  elseif maxx<1e-20; d=1e20; 
  else d=1; 
  end 
  p1=p1*d; p2=p2*d; p1s=p1s*d; p2s=p2s*d; 
end
r=p2/p2s;
