function y = PropSprod(P,Q,d,x)
% x a column n-vector, d a column n-vector, P and Q n-by-k
% y = (tril(Q*P',-1) + diag(d) + triu(P*Q',1))*x;
[n,k] = size(P);
if n<=2*k % Other reasonable nMin's acceptable
y = (tril(Q*P',-1) + diag(d) + triu(P*Q',1))*x;
else
n1 = floor(n/2);
f11 = PropSprod(P(1:n1,:),Q(1:n1,:),d(1:n1),x(1:n1));
f12 = P(1:n1,:)*(Q(n1+1:n,:)'*x(n1+1:n));
f21 = Q(n1+1:n,:)*(P(1:n1,:)'*x(1:n1));
f22 = PropSprod(P(n1+1:n,:),Q(n1+1:n,:),d(n1+1:n),x(n1+1:n));
y = [f11+f12;f21+f22];
end

% About O(k*n*logn) flops
% Full credit for correct method with this complexity.
% Any reasonable base case argument was accepted. -1 Missing discussion