function b = fastToeplitz( c, r, x ) 
% Compute b = Ax, where A is a Toeplitz matrix with its first 
% column given by the vector c and first row given by r. 
% Note r(1) should be equal to c(1). 
% 
% Author: Alex Townsend, September 2014
n = numel(x); 

if abs( r(1) - c(1) ) > eps 
    warning('Ignoring r(1).');
end

% A toeplitz matrix can be embedded into a Circulant matrix: 
% 
%  T11 = T; 
%  T12 = toeplitz([0;r(end:-1:2)],[0;c(end:-1:2)]);
%  C = [T11 T12 ; 
%       T12 T11];      % circulant
% 
%  Use this fact to compute the matrix-vector multiply using FFT. 
b = fastCirculant( [c ; 0 ; r(end:-1:2)], [x ; zeros(n,1) ] );
b = b(1:n); 

% In the 1st line, notice [v ; zeros(n,1) ], which is zeroing out the 
% contributions from the last n columns of C. 
%
% In the 2nd line, the contributions from rows n+1:2n are removed as we do not
% need them. 
% 
% Therefore, we are computing a matrix-vector product only with the principal
% nxn block (the original Toeplitz matrix). 

end