function v=ppual(pp,xx,left)
%
%      v = ppual(pp,xx[,left])
%
%  returns the value at  x  of the pp function described by  pp . The 
% function is made left-continuous if the optional third argument is present.

% copyright 1987-1993 Carl de Boor
% C. de Boor / latest change: June 19, 1989
% C. de Boor / latest change: December 1, 1990 (correct misprint in last line)
% C. de Boor / latest change: 19 se 94 (temporarily programmed around the
%                               breaks ([1 1 1 1 ... ]) discontinuity )
% C. de Boor / latest change: 19 oc 94 (optional argument for left continuity)


%  if necessary, sort  xx 
xs=xx(:)';ix=length(xs);tosort=0;
if (length(find(diff(xs)<0))>0),
   tosort=1;[xs,indx]=sort(xs);
end

%  take apart  pp
[x,c,l,k,d]=ppbrk(pp);

% for each data point, compute its breakpoint interval
if nargin<3
   index=max([sorted(x(1:l),xs);ones(1,ix)]);
else
   index=max([l+1-sorted(-x(l+1:-1:2),-xs);ones(1,ix)]);
end

% now go to local coordinates
% xs=xs-x(index);
if all(index==1), xs = xs - x(1)*index;
else xs = xs - x(index); end
% ... and apply nested multiplication:
if (d>1),
   ones(d,1)*xs; xs=ans(:)';
   1+d*ones(d,1)*index+[-d:-1]'*ones(1,ix); index=ans(:);
end
if all(index==1),
   vv=c(1,1)*index;
   for i=2:k,
      vv = xs.*vv + c(1,i)*index;
   end
else
   vv=c(index,1)';
   for i=2:k,
      vv = xs.*vv + c(index,i)';
   end
end

v=zeros(d,length(xx));v(:)=vv;
if tosort>0,[junk,indx]=sort(indx); v=v(:,indx);end