Pprod(pd)

function [Pc Pi] = Pprod(pd)
% [Pc Pi] = Pprod(pd)
% Compute multiplication of n-polynomials, and also n leave-one-out
% polynomials
%
% INPUT:
%   pd (n x m) stores n polynomials (of order m-1) in each row
% OUTPUTS:
%   - Pc is product of all polynomials: PROD [ pd(j,:) for j=1,2, ..., n ]
%   - Pi is an array polynomials of product but leaving-one-out.
%     Pi(i,:) = PROD [ pd(j,:) for j#i ]
%
% Smart booking so as the complexity is O(n*log(n))
%
% See also: StdEllPrj
%
% Author: Bruno Luong <[e-mail]>
%   Original: 24-May-2010

n = size(pd,1);
k = nextpow2(n);

% Building the tree from bottom-to-top (leaf to root)
tree = cell(1,k+1);
% initialize leaf
tree{1} = struct('child', [], 'coef', num2cell(pd,2));

k = 1;
l = n; % size(tree{k},1);
while l>1
    tk = tree{k};
    l = ceil(l/2);
    tkp1 = struct('child', cell(l,1), 'coef', cell(l,1));
    % Take convolution by pair
    for i=1:l
        j = 2*i-1;
        if j<size(tk,1)
            tkp1(i) = struct('child', [j j+1], ...
                             'coef', conv(tk(j).coef,tk(j+1).coef));
        else % alone trailing
            tkp1(i) = struct('child', j, ...
                             'coef', tk(j).coef);
        end        
    end % for-loop
    k = k+1;
    tree{k} = tkp1(:);
end % while

% Retreive the product of all the polynomials is here
Pc = tree{end}.coef;

% Building the leave-one-out products
Pi = leaveoneout(tree, cell(n, 1), size(tree,2), 1, 1);
% Put the result in array
Pi = cat(1, Pi{:});

end % Pprod

%%
% Go through the tree recursively , and take the product of polynomial
% by leaving one-out
function Pi = leaveoneout(tree, Pi, level, ind, cp)
% tree: is the tree, unchanged
% Pi: current cell list of *all* polynomials, updated when we hit
%     the first level (leaf)
% level: where we are on the tree?
% ind: indice of the node in the tree at the current level
% cp: current cumulative product of the branch currently visited

if level==1
    Pi(ind) = {cp}; % done, return
else
    % Recursive call
    node = tree{level}(ind);
    child = node.child;
    up = level-1;
    if length(child)==1
        Pi = leaveoneout(tree, Pi, up, child, cp);
    else % length(child)==1
        Pi = leaveoneout(tree, Pi, up, child(1), ...
                         conv(cp, tree{up}(child(2)).coef));
        Pi = leaveoneout(tree, Pi, up, child(2), ...
                         conv(cp, tree{up}(child(1)).coef));
    end
end

end % leaveoneout