function [p,c] = GaussChebyshevClenshaw( f_in, n_in, m_in )
% GAUSSCHEBYSHEVCLENSHAW Gauss-Chebyshev quadrature 
%                        to find approximate Chebyshev series
% GaussChebyshevClenshaw( f, n, m )
% Find a degree m Chebyshev series expansion for f on -1 <= x <= 1
% using n-poing Gauss-Chebyshev quadrature.
%
% Evaluate the resulting series using the Clenshaw algorithm at 201 points
% on -1 <= x <= 1 and compare with the original function.

if nargin<2,
    n = 40;
else
    n = n_in;
end


xi = cos( pi*((1:n)-1/2)/n );
if nargin <1,
    f = exp( - xi );
else
    f = feval( f_in, xi );
end

if nargin< 3,
    m = 15; % m+1 terms in the series we compute
else
    m = m_in+1;
end

c = zeros(m+1,1);
c(1) = 1/n *sum(f);
for k=1:m,
    c(k+1) = 2/n*f*cos( k*pi*((1:n)-1/2)/n ).';
end

t = linspace(-1,1,201);

function p = ClenshawChebyshev( t, c )
% Clenshaw algorithm for evaluation of the series
m = length(c) - 1;
if m<0,
    % empty c ==> zero polynomial
    p = zeros(size(t));
elseif m==0,
    % degree 0 ==> constant polynomial
    p = c(1)*ones(size(t));
else
    % degree 1 or more, use Clenshaw algorithm
    y2 = zeros(size(t));  % y_{m+1}
    y1 = c(m+1)*ones(size(t));  % y_{m}
    % 
    % p = sum( c(j)*T(j-1,t), j=1..m-1) + c(m)*T(m-1,t) 
    %        + y1*T(m,t) - y2*T(m-1,t)  Loop invariant initialization
    for k=m:-1:2,
        % p = sum( c(j)*T(j-1,t), j=1..k-1) + c(k)*T(k-1,t) + y1*T(k,t) 
        %  - y2*T(k-1,t)                Loop invariant
        y0 = c(k) +2*t.*y1 - y2;   % y_{k-1} 
        % p = sum( c(j)*T(j-1,t), j=1..k-1) + y0*T(k-1,t) +
        % y1*(T(k,t)-2tT(k-1,t)) [= -y1*T(k-2,t) ]
        y2 = y1; 
        y1 = y0; 
        % p = sum( c(j)*T(j-1,t), j=1..k-2) + c(k-1)*T(k-2,t) + 
        % y1*T(k-1,t) - y2*T(k-2,t)     Proving correctness
    end
    % p = c(1)*T(0,t) + y1*T(1,t) - y2*T(0,t)
    p = c(1)-y2 + y1.*t;
end
end

p = ClenshawChebyshev( t, c );
if nargin<1,
    ft = exp(-t);
    errs = p - ft ;
else
    ft = feval(f_in,t);
    errs = p - ft;
end
figure(1), plot( t, errs, 'k'), set(gca,'fontsize',16),xlabel('t'),ylabel('absolute error')
figure(2), plot( t, p, 'k', t, ft, 'k--'), set(gca,'fontsize',16),xlabel('t'),ylabel('Chebyshev series')
end