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]);