function v = mpval(p,x)
%MPVAL evaluate multivariate polynomial in (shifted, normalized) power form
%
%        V = MPVAL(P,X)
%
% Returns the value at  X  of the (d-variate) polynomial whose (shifted,
% normalized) power form 
%
%   p(x) = sum_a (x-center)^a (|a| choose a) coefs(:,nn(a))
%
% is contained in  P , with  X  of size  d-by-nx .
%
% The powers are so normalized that a multivariate version of Horner's scheme
% (aka the de Casteljau algorithm) works simply; see below.
%
% See MPMAK for details on the form.

% cb: 26jul97, 9aug97; reverted to v4 (horrible) 7-13sep97
%   Copyright (c) 1997 by C. de Boor

[dx,nx] = size(x);
[coefs,d,k,center,dimpjd,mm] = mpbrk(p);
[df,nc] = size(coefs);

if dx~=d 
   error(['The given polynomial is ',int2str(d),'-variate', ...
            ', while  X  is in R^',int2str(dx),'.'])
end

x = x - center(:,ones(1,nx));

%   Horner's scheme (aka Nested Multiplication) consists of the following loop:
%
%     vk <-- (coefs(nn(a)) : |a|=k)
%     for j=k-1:-1:0
%        vj <-- (coefs(nn(a)) : |a|=j) + sum_{i=1:d} x(i)*vjp1(n(a+e_i))
%     end
%
%   Hence  v0 = sum_a coefs(nn(a)) x^a * # sequences (0=b_0,b_1,...,b_|a|=a) 
%   with  b_{j+1}-b_j in {e_1,...,e_d}, all j . In other words,
%       v0 = sum_a x^a (|a| choose a) coefs(nn(a)) = p(x) .
%
%   For this, it is important to know that
%   DIMPJD(2+j) = dim Pi_j(R^d) = ( j + d choose d ) is the number of d-variate 
%      monomials of degree  <=j , while   
%   M(i,nn(a)) = n(a+e_i), i=1:d , and
%   nn(a) = dim Pi_{<|a|} + n(a) .

dd = [1:d].'; v = zeros(df,nx);
for ff=1:df
   coefs(ff,dimpjd(k+1)+1:dimpjd(k+2)).'; vff = ans(:,ones(1,nx));
   if k>0
      for j=k-1:-1:0
         rangej = dimpjd(j+1)+1:dimpjd(j+2); lj = length(rangej);
         reshape(coefs(ff,rangej),lj,1);
         if d>1
            vff = ans(:,ones(1,nx)) + ...
                reshape(sum(reshape(vff(mm(:,rangej),:),d,lj*nx).* ...
                        reshape(x(dd(:,ones(1,lj)),:),d,lj*nx)),lj,nx);
            else
            vff = ans(:,ones(1,nx)) + ...
                vff(mm(:,rangej),:).*x(dd(:,ones(1,lj)),:);
         end
      end
   end
   v(ff,:) = reshape(vff,1,nx);
end

%v5 v = coefs(:,dimpjd(k+1)+1:dimpjd(k+2),ones(1,nx));
%v5 if k>0
%v5    x = reshape(x,[1,d,1,nx]);
%v5    for j=k-1:-1:0
%v5       rangej = dimpjd(j+1)+1:dimpjd(j+2); lj = length(rangej);
%v5       v = coefs(:,rangej,ones(1,nx)) + ...
%v5               reshape(sum(reshape(v(:,mm(:,rangej),:),[df,d,lj,nx]).* ...
%v5                   x(ones(df,1),:,ones(lj,1),:),2),[df,lj,nx]);
%v5    end
%v5 end
%v5 v = reshape(v,[df,nx]);