--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,12 @@
+2009-11-02  Jaroslav Hajek  <highegg@gmail.com>
+
+	* polynomial/mkpp.m: Improve Matlab compatibility.
+	* polynomial/ppval.m: Ditto.
+	* polynomial/unmkpp: Update.
+	* polynomial/pchip.m: Update and optimize.
+	* polynomial/spline.m: Update.
+	* general/__splinen__.m: Update.
+
 2009-10-23  Jaroslav Hajek  <highegg@gmail.com>
 
 	* general/tril.m, general/triu.m: Remove sources.

--- a/scripts/general/__splinen__.m
+++ b/scripts/general/__splinen__.m
@@ -37,7 +37,7 @@
   endif
   yi = y;
   for i = length(x):-1:1
-    yi = permute (spline (x{i}, yi, xi{i}), [length(x),1:length(x)-1]);
+    yi = permute (spline (x{i}, yi, xi{i}(:)), [length(x),1:length(x)-1]);
   endfor
 
   [xi{:}] = ndgrid (cellfun (@(x) x(:), xi, "UniformOutput", false){:});

--- a/scripts/polynomial/mkpp.m
+++ b/scripts/polynomial/mkpp.m
@@ -27,14 +27,10 @@
 ## lowest.  There must be one row for each interval in @var{x}, so 
 ## @code{rows (@var{p}) == length (@var{x}) - 1}.  
 ##
-## You can concatenate multiple polynomials of the same order over the 
-## same set of intervals using @code{@var{p} = [ @var{p1}; @var{p2}; 
-## @dots{}; @var{pd} ]}.  In this case, @code{rows (@var{p}) == @var{d} 
-## * (length (@var{x}) - 1)}.
-##
-## @var{d} specifies the shape of the matrix @var{p} for all except the
-## last dimension.  If @var{d} is not specified it will be computed as
-## @code{round (rows (@var{p}) / (length (@var{x}) - 1))} instead.
+## @var{p} may also be a multi-dimensional array, specifying a vector-valued
+## or array-valued polynomial. The shape is determined by @var{d}. If @var{d} is
+## not given, the default is @code{size (p)(1:end-2)}. If @var{d} is given, the
+## leading dimensions of @var{p} are reshaped to conform to @var{d}.
 ##
 ## @seealso{unmkpp, ppval, spline}
 ## @end deftypefn
@@ -44,14 +40,24 @@
     print_usage ();
   endif
   pp.x = x(:);
-  pp.P = P;
-  pp.n = length (x) - 1;
-  pp.k = columns (P);
+  n = length (x) - 1;
+  if (n < 1)
+    error ("mkpp: at least one interval is needed");
+  endif
+  nd = ndims (P);
+  k = size (P, nd);
   if (nargin < 3)
-    d = round (rows (P) / pp.n); 
+    if (nd == 2)
+      d = 1;
+    else
+      d = prod (size (P)(1:nd-1));
+    endif
   endif
   pp.d = d;
-  if (pp.n * prod (d) != rows (P))
+  pp.P = P = reshape (P, prod (d), [], k);
+  pp.orient = 0;
+
+  if (size (P, 2) != n)
     error ("mkpp: num intervals in x doesn't match num polynomials in P");
   endif
 endfunction

--- a/scripts/polynomial/pchip.m
+++ b/scripts/polynomial/pchip.m
@@ -115,7 +115,7 @@
   c2 = -c3 - del1;
   c3 = c3 / h;
 
-  coeffs = [c3.'(:), c2.'(:), d1.'(:), f1.'(:)];
+  coeffs = cat (3, c3, c2, d1, f1);
   pp = mkpp (x, coeffs, szy(1:end-1));
 
   if (nargin == 2)

--- a/scripts/polynomial/ppval.m
+++ b/scripts/polynomial/ppval.m
@@ -19,9 +19,11 @@
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi})
 ## Evaluate piece-wise polynomial @var{pp} at the points @var{xi}.  
-## If @code{@var{pp}.d} is a scalar greater than 1, or an array,
-## then the returned value @var{yi} will be an array that is 
-## @code{d1, d1, @dots{}, dk, length (@var{xi})]}.
+## If @var{pp} is scalar-valued, the result is an array of the same shape as @var{xi}.
+## Otherwise, the size of the result is @code{[pp.d, length(@var{xi})]} if
+## @var{xi} is a vector, or @code{[pp.d, size(@var{xi})]} if it is a multi-dimensional
+## array. If pp.orient is 1, the dimensions are permuted as in interp1, to
+## @code{[pp.d, length(@var{xi})]} and @code{[pp.d, size(@var{xi})]} respectively.
 ## @seealso{mkpp, unmkpp, spline}
 ## @end deftypefn 
 
@@ -33,23 +35,45 @@
   if (! isstruct (pp))
     error ("ppval: expects a pp structure");
   endif
-  if (isempty (xi))
-    yi = [];
-  else
-    transposed = (columns (xi) == 1);
-    xi = xi(:);
-    xn = length (xi);
-    idx = lookup (pp.x, xi, "lr");
-    dx = (xi - pp.x(idx)).';
-    dx = reshape (dx(ones(1,prod(pp.d)),:),[pp.d,xn]);
-    c = reshape (pp.P(:,1), pp.n, prod (pp.d));
-    yi = reshape (c(idx,:).', [pp.d, xn]);
-    for i  = 2 : pp.k;
-      c = reshape (pp.P(:,i), pp.n, prod (pp.d));
-      yi = yi .* dx + reshape (c(idx,:).', [pp.d, xn]);
-    endfor
-    if (transposed && isscalar (pp.d) && pp.d == 1)
-      yi = yi.';
-    endif
+
+  ## Extract info.
+  x = pp.x;
+  P = pp.P;
+  d = pp.d;
+  k = size (P, 3);
+  nd = size (P, 1);
+
+  ## Determine resulting shape.
+  if (d == 1) # scalar case
+    yisz = size (xi);
+  elseif (isvector (xi)) # this is special
+    yisz = [d, length(xi)];
+  else # general
+    yisz = [d, size(xi)];
   endif
+
+  ## Determine intervals.
+  xi = xi(:);
+  xn = numel (xi);
+
+  idx = lookup (x, xi, "lr");
+
+  ## Offsets.
+  dx = (xi - x(idx)).';
+  dx = dx(ones (1, nd), :); # spread (do nothing in 1D)
+
+  ## Use Horner scheme.
+  yi = P(:,idx,1);
+  for i = 2:k;
+    yi .*= dx;
+    yi += P(:,idx,i);
+  endfor
+
+  ## Adjust shape.
+  yi = reshape (yi, yisz);
+  if (d != 1 && pp.orient == 1)
+    ## Switch dimensions to match interp1 order.
+    yi = shiftdim (yi, length (d));
+  endif
+
 endfunction

--- a/scripts/polynomial/spline.m
+++ b/scripts/polynomial/spline.m
@@ -186,7 +186,7 @@
   c = c(1:n-1,:);
   b = b(1:n-1,:);
   a = a(1:n-1,:);
-  coeffs = [d(:), c(:), b(:), a(:)];
+  coeffs = cat (3, d.', c.', b.', a.');
   ret = mkpp (x, coeffs, szy(1:end-1));
 
   if (nargin == 3)

--- a/scripts/polynomial/unmkpp.m
+++ b/scripts/polynomial/unmkpp.m
@@ -30,10 +30,7 @@
 ## (@var{i}, :)} contains the coefficients for the polynomial over
 ## interval @var{i} ordered from highest to lowest.  If @code{@var{d} >
 ## 1}, @code{@var{p} (@var{r}, @var{i}, :)} contains the coefficients for 
-## the r-th polynomial defined on interval @var{i}.  However, this is 
-## stored as a 2-D array such that @code{@var{c} = reshape (@var{p} (:,
-## @var{j}), @var{n}, @var{d})} gives @code{@var{c} (@var{i},  @var{r})}
-## is the j-th coefficient of the r-th polynomial over the i-th interval.
+## the r-th polynomial defined on interval @var{i}.
 ## @item @var{n}
 ## Number of polynomial pieces.
 ## @item @var{k}
@@ -54,7 +51,10 @@
   endif
   x = pp.x;
   P = pp.P;
-  n = pp.n;
-  k = pp.k;
+  n = size (P, 2);
+  k = size (P, 3);
   d = pp.d;
+  if (d == 1)
+    P = reshape (P, n, k);
+  endif
 endfunction