function v=mpolyval(p,x)
% v=mpolyval(p,x)
% evaluate multivariate polynomial
% v =sum p(i_1,...,i_m)*x(1)^(n_1-i_1)...x(m)^(n_m-i_m)
% using recursively the Horners rule
% note: if p represents a univariate polynomial p must be a column

n=size(p);
m=length(n);
% delete trailing 1 if 2d matrix
if m==2
    if n(2)==1,
        n(2)=[]; 
        m=1;
    end
end

if m==1,            % if 1d, evaluate directly
    v=polyval(p,x(1)); 
else                % recursive
    q=zeros(n(m),1);
    p2=reshape(p,[prod(n(1:m-1)),n(m)]);
    for i=1:n(m),
        psel=reshape(p2(:,i),[n(1:m-1),1]);
        q(i)=mpolyval(psel,x(1:m-1));
    end
    v=polyval(q,x(m));
end
%disp('mpolyval:'),p,x,v