  function P1()
% Checks out FastProd

% Small example to check correctness
t = 10; n = 2^t; d = randn(n,1); S = randn(n,2); T = randn(n,2);
A =  tril(S*T',-1) + diag(d) + triu(T*S',1);
x = randn(n,1);
y0 = A*x;
y = FastProd(d,S,T,x);
clc
fprintf('n=%1d\n  ||y - Ax||/||Ax|| = %10.3e\n',n,norm(y-y0)/norm(y0))

% Large example to check efficiecy
t = 15; n = 2^t; d = randn(n,1); S = randn(n,2); T = randn(n,2); x = randn(n,1);
tic
yCVL0 = FastProdCVL0(d,S,T,x);
tCVL0 = toc;
tic
yCVL1 = FastProdCVL1(d,S,T,x);
tCVL1 = toc;
tic
yCVL2 = FastProdCVL2(d,S,T,x);
tCVL2 = toc;
tic
yYOU = FastProd(d,S,T,x);
tYOU = toc;
relErr = norm(yYOU-yCVL0)/norm(yCVL0);
fprintf('\nn=%1d, tYOU/tCVL0 = %6.1f, relErr = %10.3e\n',n,tYOU/tCVL0,relErr)
fprintf('\nn=%1d, tYOU/tCVL1 = %6.1f, relErr = %10.3e\n',n,tYOU/tCVL1,relErr)
fprintf('\nn=%1d, tYOU/tCVL2 = %6.1f, relErr = %10.3e\n',n,tYOU/tCVL2,relErr)




% This O(n) implementation is based on spotting repetitive pointwise vector
% products.

   function y = FastProdCVL0(d,S,T,x)
% d and x are column n vectors, S and T are n-by-2, n is a positive power of two.
% y = A*x where A = diag(d) + tril(S*T',-1) + triu(T*S',1)
n = length(d);
y = d.*x + ...
   [0;S(2:n,1).*cumsum(T(1:n-1,1).*x(1:n-1))] + ...
   [0;S(2:n,2).*cumsum(T(1:n-1,2).*x(1:n-1))] + ...
   [T(1:n-1,1).*flip(cumsum(flip(S(2:n,1).*x(2:n))));0] + ...
   [T(1:n-1,2).*flip(cumsum(flip(S(2:n,2).*x(2:n))));0];

% This O(nlogn) recursive implementation is based on divide and conquer..

   function y = FastProdCVL1(d,S,T,x)
% d and x are column n vectors, S and T are n-by-2, n is a positive power of two.
% y = A*x where A = diag(d) + tril(S*T',-1) + triu(T*S',1)
 n = length(d);
 if n==256
     % Base case is done conventionally...
     A0 = tril(S*T',-1);
     y = (A0+A0')*x + d.*x;
 else
    m = n/2;
 %  Partition:
 %       S1 = S(1:m,:)     S2 = S(m+1:n,:)
 %       T1 = T(1:m,:)     T2 = T(m+1:n,:)
 %       D1 = diag(d(1:m)) D2 = diag(d(m+1:n))
 %       x1 = x(1:m)       x2 = x(m+1:n)
 %       y1 = y(1:m)       y2 = y(m+1:n)
 %       A1 = tril(S1*T1',-1)+D1+triu(T1*S1',1)
 %       A2 = tril(S2*T2',-1)+D2+triu(T2*S2',1)
 %       
 %  [y1]     [   A1     T1*S2' ][x1]   [ A1*x1 + T1*(S2'*x2) ]
 %  [  ]   = [                 ][  ] = [                     ]
 %  [y2]     [ S2*T1'    A2    ][x2]   [ S2*(T1'*x1) + A2*x2 ]           
     m = n/2;
     y1 = FastProdCVL1(d(1:m),S(1:m,:),T(1:m,:),x(1:m)) + ...
                                        T(1:m,:)*(S(m+1:n,:)'*x(m+1:n));
     y2 = FastProdCVL1(d(m+1:n),S(m+1:n,:),T(m+1:n,:),x(m+1:n)) +...
                                        S(m+1:n,:)*(T(1:m,:)'*x(1:m));
     y = [y1;y2];
 end% 
 
% This O(nlogn) implementation is a nonrecursive implementation of
% FastProdCL0

  function y = FastProdCVL2(d,S,T,x)
% d and x are column n vectors, S and T are n-by-2, n is a positive power of two.
% y = A*x where A = diag(d) + tril(S*T',-1) + triu(T*S',1)
n = length(d);
t = log2(n);
y = d.*x;
for k=1:t
   jdx = 2^(k-1)+1:2^k;
   idx = jdx-2^(k-1);
   for j=1:n/2^k
       y(idx) = y(idx)+T(idx,:)*(S(jdx,:)'*x(jdx));
       y(jdx) = y(jdx)+S(jdx,:)*(T(idx,:)'*x(idx));
       jdx = jdx + 2^k;
       idx = jdx - 2^(k-1);
   end
end
