function polvalues=polynomials(points,order)
% generates all polynomials on points up to order
% The order MUST be increasing from left to right.
[numpoints,dim]=size(points); % handle trivial cases first
if order==0
    polvalues=[];
    return
end
if order==1
    polvalues=ones(numpoints,1);
    return
end
if order==2
    polvalues=[ones(numpoints,1) points];
    return
end       
if dim==1
  polvalues=zeros(numpoints,order);
  polvalues(:,1)=ones(numpoints,1);
  polvalues(:,2)=points(:,1);
  for i=3:order
      polvalues(:,i)=polvalues(:,i-1).*points(:,1);
  end
  return
end           % general case done recursively
polvalues=[polynomials(points,order-1)...
        polynomials_exact_order(points,order)];
