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