function [ froots, Einter, root_residuals ] = cpr ( f, a, b, N )

%*****************************************************************************80
%
%% cpr() computes the real roots of f(x) on the real interval [a,b].
%
%  Discussion:
%
%    A Chebyshev Proxy Rootfinder is used.
%
%    f(x) is assumed to be a smooth real-valued scalar function of a scalar
%    argument.
%
%    The code will only search for real roots of f(x), within the
%    user-specified interval.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 March 2024
%
%  Author:
%
%    John Boyd
%
%  Reference:
%
%    John Boyd,
%    Solving Transcendental Equations,
%    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
%    Perturbation Series, and Oracles,
%    SIAM, 2014,
%    ISBN: 978-1-611973-51-8,
%    LC: QA:353.T7B69
%
%  Input:
%
%    string f: the name of the function.
%
%    real a, b: the endpoints of the search interval.
%
%    integer N: the degree of Chebyshev expansion, using N+1 points.
%
%  Output:
%
%    real froots(nroots): a column vector, contains the computed roots.
%
%    real Einter: ???
%
%    real root_residuals(nroots): a column vector, the function value at 
%    each computed root.
%
%  Local:
%
%    real dyn_range_tol: trigger for warning of small dynamic range.
%
%    real epscutoff: trailing coefficient j is discarded if 
%    abs ( acoeff(j)) < epscutoff * max ( abs ( acoeff() ) ).
%
%    real sigma: Discard any root wih ( 1 + sigma ) < abs ( Re ( root ) )
%    measured in Chebyshev coordinate range [-1,+1].
%
%    real tau: Discard any root with tau < abs ( Im ( root ) ).
%
  epscutoff = 1.0E-13;
  tau = 1.0E-08;
  sigma = 1.0E-06;
  dyn_range_tol = 1.0E-08;
%
%  Compute trigonometic argument row vector T.
%
  t = ( 0 : 1 : N ) * pi / N;
% t = linspace ( 0.0, pi, N + 1 );
%
%  Set row vector of Chebyshev nodes in [-1,+1].
%
  xi = cos ( t );
%
%  Compute row vector of interpolation nodes in [a,b].
%
  x = 0.5 * ( b - a ) * xi + ( b + a ) * 0.5;
%
%  Create column vector of function samples.
%  Since f might not be vectorizable, we use a loop.
%
  fa = zeros ( 1, N + 1 );
  for i = 1 : N + 1
    fa(i) = f ( x(i) );
  end
%
%  Compute acoeff, the Chebyshev coefficients.
%
  phia = zeros ( N + 1, N + 1 );
  pj = ones ( 1, N + 1 );
  pj(1) = 2.0;
  pj(N+1) = 2.0;
  for i = 1 : N + 1
    phia(i,:) = ( 2.0 / N / pj(i) ) * cos ( ( i - 1 ) * t(:) ) ./ pj(:);
  end
  acoeff = phia * fa';
%
%  Truncate the tail of the series, if very small.
%  
  acoeffmax = max ( abs ( acoeff ) );
  Nt = N;
  tailnorm = 0.0;
  for k = N : -1 : 2
    tailnorm = tailnorm + abs ( acoeff(k+1) );
    if ( tailnorm < epscutoff * acoeffmax )
      Nt = Nt - 1;
    end
  end
%
%  Create the Chebyshev companion matrix.
%
  A = zeros ( Nt, Nt );
  A(1,2) = 1.0;
  for j = 2 : Nt - 1
    A(j,j-1) = 0.5;
    A(j,j+1) = 0.5;
  end
  A(Nt,1:Nt) = - acoeff(1:Nt) / ( 2.0 * acoeff(Nt+1) );
  A(Nt,Nt-1) = A(Nt,Nt-1) + 0.5;
%
%  eig() returns a column vector of eigenvalues.
%  The roots are in Chebyshev coordinates xi.
%
  all_roots_in_xi = eig ( A );
%
%  Accept roots only if within tau of real axis and Chebyshev interval.
%
  nroots = 0;
  for i = 1 : Nt
    evi = all_roots_in_xi(i);
    if ( abs ( imag ( evi ) ) < tau * abs ( real ( evi ) ) )
      if ( abs ( real ( evi ) ) <= 1.0 + sigma )
        nroots = nroots + 1;
        froots(nroots,1) = 0.5 * ( b - a ) * real ( evi ) + ( b + a ) / 2.0;
      end
    end
  end
%
%  Sort the roots.
%
  froots = sort ( froots );
%
%  Compute column vector of function residuals at computed roots.
%
  root_residuals = f ( froots );
%
%  Interstitial interpolation residual.
%
  tinter = t(1:N) + 0.5 / N;
  xiint = cos ( tinter );
  xinter = 0.5 * ( b - a ) * xiint + ( b + a ) * 0.5;
  fainter = zeros ( 1, N );
  for i = 1 : N
    fainter(i) = f ( xinter(i) );
  end

  taall = [ t, tinter ];
  faall = [ fa, fainter ];
  xiall = [ xi, xiint ];
  xall = [ x, xinter ];

  fpoly = zeros ( 1, 2 * N + 1 );
  for j = 1 : N + 1
    fpoly = fpoly + acoeff(j) * cos ( ( j - 1 ) * taall );
  end
  Einter = max ( abs ( faall((N+2):(2*N+1)) - fpoly((N+2):(2*N+1)) ) );
%
%  Compute dynamic range.
%
  fmax = max ( abs ( faall ) );
  dyn_range = min ( abs ( faall ) ) / fmax;
  if ( dyn_range < dyn_range_tol )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'cpr() - Warning!\n' );
    fprintf ( 1, '  dynamic range = %g < %g\n', dyn_range, dyn_range_tol );
  end

  return
end

