Mercurial > hg > octave-lojdl

changeset 12623:b0e60ad4ae26

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge stable to default
author John W. Eaton <jwe@octave.org>
date Tue, 19 Apr 2011 21:23:28 -0400
parents 32d5186266ab (diff) 6f3f18957851 (current diff)
children b1f4f0eba971
files
diffstat 43 files changed, 698 insertions(+), 535 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog
+++ b/ChangeLog
@@ -2,6 +2,10 @@
 
 	* NEWS: Add colstyle to list of new functions for 3.4
 
+2011-04-08  Rik  <octave@nomad.inbox5.com>
+
+	* NEWS: Deprecate studentize(), add new function zscore().
+
 2011-04-04  Rik  <octave@nomad.inbox5.com>
 
 	* NEWS: Add perror, strerror to list of functions deprecated in 3.4
--- a/NEWS
+++ b/NEWS
@@ -7,7 +7,8 @@
 
     iscolumn
     issrow
- 
+    zscore
+
  ** Deprecated functions.
 
     The following functions were deprecated in Octave 3.2 and have been
@@ -30,7 +31,8 @@
     release after 3.6):
 
       is_duplicate_entry
-      
+      studentize
+
 
 Summary of important user-visible changes for version 3.4:
 ---------------------------------------------------------
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,7 @@
+2011-04-16  Ben Abbott  <bpabbott@mac.com>
+
+	* interpreter/plot.txi: Clarify that inheritance of default property
+	values only applies to the named object type.
 2011-04-14  Rik  <octave@nomad.inbox5.com>
 
 	* interpreter/plot.txi: Add colstyle function to documentation.
@@ -16,6 +20,11 @@
 	* interpreter/func.txi: Add discussion of isargout to Ignoring
 	Arguments section of documentation.
 
+2011-04-08  Rik  <octave@nomad.inbox5.com>
+
+	* interpreter/stats.txi: Deprecate studentize(), replace with
+	zscore().
+
 2011-04-07  Rik  <octave@nomad.inbox5.com>
 
 	* interpreter/Makefile.am: Add spellcheck target to documentation
@@ -42,7 +51,7 @@
 
 2011-04-04  Rik  <octave@nomad.inbox5.com>
 
-	* interpreter/doccheck/aspell-octave.en.pws, interpreter/nonlin.txi, 
+	* interpreter/doccheck/aspell-octave.en.pws, interpreter/nonlin.txi,
 	interpreter/tips.txi: Spellcheck documentation for 3.4.1 release.
 
 2011-04-04  Rik  <octave@nomad.inbox5.com>
--- a/doc/interpreter/plot.txi
+++ b/doc/interpreter/plot.txi
@@ -2390,7 +2390,11 @@
 may override the factory defaults.
 
 Although default values may be set for any object, they are set in
-parent objects and apply to child objects.  For example,
+parent objects and apply to child objects, of the specified object type.
+For example, seeting the default @code{color} property of @code{line}
+objects to "green", for the @code{root} object, will result in all
+@code{line} objects inheriting the @code{color} "green" as the default
+value.
 
 @example
 set (0, "defaultlinecolor", "green");
--- a/doc/interpreter/stats.txi
+++ b/doc/interpreter/stats.txi
@@ -114,7 +114,7 @@
 
 @DOCSTRING(center)
 
-@DOCSTRING(studentize)
+@DOCSTRING(zscore)
 
 @DOCSTRING(histc)
 
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,7 @@
+2011-04-12  Rik  <octave@nomad.inbox5.com>
+
+	* LSODE.cc: Add semicolon to error messages to prevent run-together text.
+
 2011-04-01  Jordi Gutiérrez Hermoso  <jordigh@gmail.com>
 
 	* MatrixType (MatrixType::operator =): Plug memory leak due to
--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -292,7 +292,7 @@
     case -2:  // excess accuracy requested (tolerances too small).
     case -3:  // invalid input detected (see printed message).
     case -4:  // repeated error test failures (check all inputs).
-    case -5:  // repeated convergence failures (perhaps bad jacobian
+    case -5:  // repeated convergence failures (perhaps bad Jacobian
               // supplied or wrong choice of mf or tolerances).
     case -6:  // error weight became zero during problem. (solution
               // component i vanished, and atol or atol(i) = 0.)
@@ -349,13 +349,13 @@
 
     case -4:
       retval = std::string ("repeated error test failures (t = ")
-        + t_curr + "check all inputs)";
+        + t_curr + "; check all inputs)";
       break;
 
     case -5:
       retval = std::string ("repeated convergence failures (t = ")
         + t_curr
-        + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)";
+        + "; perhaps bad Jacobian supplied or wrong choice of integration method or tolerances)";
       break;
 
     case -6:
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,8 +1,47 @@
+2011-04-18  Paul Boven  <p.boven@xs4all.nl>
+
+	* image/image.m: Fixed naming of variables in texinfo
+
+2011-04-17  Patrick Häcker  <magicmuscleman>
+
+	* strings/mat2str.m: Limit the number of digits to one less than
+	available for double.
+
+2011-04-15  Kai Habel  <kai.habel@gmx.de>
+
+	* general/interp1.m, polynomial/mkpp.m, polynomial/pchip.m,
+	polynomial/ppder.m, polynomial/ppint.m, polynomial/ppjumps.m,
+	polynomial/ppval.m, polynomial/spline.m, polynomial/unmkpp.m:
+	Make functions more compatible with respect to handling of
+	picewise polynoms (pp). Rename pp-struct elements.
+	Handle nD-arguments correctly. Tests added.
+	(bugs #32040, #32045)
+
 2011-04-13  David Bateman  <dbateman@free.fr>
 
 	* plot/colstyle.m : New function.
 	* plot/module.mk plot_FCN_FILES) : Add it here.
 
+2011-04-13  Rik  <octave@nomad.inbox5.com>
+
+	* help/__makeinfo__.m: Simplify function by using regular expressions.
+	Eliminate third input argument see_also function.
+
+2011-04-13  Rik  <octave@nomad.inbox5.com>
+
+	* general/isdir.m, general/isequal.m, general/isequalwithequalnans.m,
+	general/isscalar.m, general/issquare.m, general/isvector.m: Refactor
+	code to put input validation first.
+
+	* general/iscolumn.m, general/isrow.m : Remove useless initialization
+	of output variable.
+
+	* general/isa.m: Add additional tests for various classes.
+
+2011-04-13  Rik  <octave@nomad.inbox5.com>
+
+	* ChangeLog: Remove results of bad merge in ChangeLog.
+
 2011-04-12  Ben Abbott  <bpabbott@mac.com>
 
 	* miscellaneous/getappdata.m: If no property name is provided, return
@@ -20,10 +59,21 @@
 
 2011-04-08  Rik  <octave@nomad.inbox5.com>
 
+	* deprecated/module.mk, statistics/base/center.m,
+	statistics/base/module.mk: Deprecate studentize(), replace with
+	zscore().
+
+2011-04-08  Rik  <octave@nomad.inbox5.com>
+
 	* linear-algebra/cond.m, linear-algebra/expm.m, linear-algebra/logm.m,
 	linear-algebra/null.m, linear-algebra/orth.m, linear-algebra/rank.m,
 	linear-algebra/rref.m: Improve docstrings.
 
+2011-04-08  Rik  <octave@nomad.inbox5.com>
+
+	* statistics/base/mode.m, statistics/base/quantile.m: Return output
+	of same class as input.
+
 2011-04-06  Rik  <octave@nomad.inbox5.com>
 
 	* miscellaneous/pack.m: Improve docstring.
--- a/scripts/deprecated/module.mk
+++ b/scripts/deprecated/module.mk
@@ -23,6 +23,7 @@
   deprecated/sphcat.m \
   deprecated/spvcat.m \
   deprecated/strerror.m \
+  deprecated/studentize.m \
   deprecated/values.m \
   deprecated/weibcdf.m \
   deprecated/weibinv.m \
copy from scripts/statistics/base/studentize.m
copy to scripts/deprecated/studentize.m
--- a/scripts/statistics/base/studentize.m
+++ b/scripts/deprecated/studentize.m
@@ -32,6 +32,12 @@
 ## Description: Subtract mean and divide by standard deviation
 
 function t = studentize (x, dim)
+  persistent warned = false;
+  if (! warned)
+    warned = true;
+    warning ("Octave:deprecated-function",
+             "studentize is obsolete and will be removed from a future version of Octave; please use zscore instead");
+  endif
 
   if (nargin != 1 && nargin != 2)
     print_usage ();
--- a/scripts/general/dblquad.m
+++ b/scripts/general/dblquad.m
@@ -37,22 +37,23 @@
 ##
 ## The optional argument @var{quadf} specifies which underlying integrator
 ## function to use.  Any choice but @code{quad} is available and the default
-## is @code{quadgk}.
+## is @code{quadcc}.
 ##
 ## Additional arguments, are passed directly to @var{f}.  To use the default
-## value for @var{tol} or @var{quadf} one may pass an empty matrix ([]).
+## value for @var{tol} or @var{quadf} one may pass ':' or an empty matrix ([]).
 ## @seealso{triplequad, quad, quadv, quadl, quadgk, quadcc, trapz}
 ## @end deftypefn
 
-function q = dblquad(f, xa, xb, ya, yb, tol, quadf, varargin)
+function q = dblquad (f, xa, xb, ya, yb, tol = 1e-6, quadf = @quadcc, varargin)
+
   if (nargin < 5)
     print_usage ();
   endif
-  if (nargin < 6 || isempty (tol))
+  if (isempty (tol))
     tol = 1e-6;
   endif
-  if (nargin < 7 || isempty (quadf))
-    quadf = @quadgk;
+  if (isempty (quadf))
+    quadf = @quadcc;
   endif
 
   inner = @__dblquad_inner__;
@@ -72,10 +73,10 @@
   endfor
 endfunction
 
-%% Nasty integrand to show quadgk off
+%% Nasty integrand to show quadcc off
 %!assert (dblquad (@(x,y) 1 ./ (x+y), 0, 1, 0, 1), 2*log(2), 1e-6)
 
-%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [],  @quadgk), pi * erf(1).^2, 1e-6)
-%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [],  @quadl), pi * erf(1).^2, 1e-6)
-%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, [],  @quadv), pi * erf(1).^2, 1e-6)
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6, @quadgk), pi * erf(1).^2, 1e-6)
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6, @quadl), pi * erf(1).^2, 1e-6)
+%!assert (dblquad (@(x,y) exp(-x.^2 - y.^2) , -1, 1, -1, 1, 1e-6, @quadv), pi * erf(1).^2, 1e-6)
 
--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -43,7 +43,7 @@
 ## Piecewise cubic Hermite interpolating polynomial
 ##
 ## @item 'cubic'
-## Cubic interpolation from four nearest neighbors
+## Cubic interpolation (same as @code{pchip})
 ##
 ## @item 'spline'
 ## Cubic spline interpolation---smooth first and second derivatives
@@ -112,7 +112,7 @@
   method = "linear";
   extrap = NA;
   xi = [];
-  pp = false;
+  ispp = false;
   firstnumeric = true;
 
   if (nargin > 2)
@@ -123,7 +123,7 @@
         if (strcmp ("extrap", arg))
           extrap = "extrap";
         elseif (strcmp ("pp", arg))
-          pp = true;
+          ispp = true;
         else
           method = arg;
         endif
@@ -138,7 +138,7 @@
     endfor
   endif
 
-  if (isempty (xi) && firstnumeric && ! pp)
+  if (isempty (xi) && firstnumeric && ! ispp)
     xi = y;
     y = x;
     x = 1:numel(y);
@@ -150,9 +150,8 @@
   szx = size (xi);
   if (isvector (y))
     y = y(:);
-  elseif (isvector (xi))
-    szx = length (xi);
   endif
+  
   szy = size (y);
   y = y(:,:);
   [ny, nc] = size (y);
@@ -191,147 +190,85 @@
 
   switch (method)
   case "nearest"
-    if (pp)
-      yi = mkpp ([x(1); (x(1:nx-1)+x(2:nx))/2; x(nx)], y, szy(2:end));
+    pp = mkpp ([x(1); (x(1:nx-1)+x(2:nx))/2; x(nx)], shiftdim (y, 1), szy(2:end));
+    pp.orient = "first";
+    
+    if (ispp)
+      yi = pp;
     else
-      idx = lookup (0.5*(x(1:nx-1)+x(2:nx)), xi) + 1;
-      yi = y(idx,:);
+      yi = ppval (pp, reshape (xi, szx));
     endif
   case "*nearest"
-    if (pp)
-      yi = mkpp ([x(1); x(1)+[0.5:(nx-1)]'*dx; x(nx)], y, szy(2:end));
+    pp = mkpp ([x(1), x(1)+[0.5:(nx-1)]*dx, x(nx)], shiftdim (y, 1), szy(2:end));
+    pp.orient = "first";
+    if (ispp)
+      yi = pp;
     else
-      idx = max (1, min (ny, floor((xi-x(1))/dx+1.5)));
-      yi = y(idx,:);
+      yi = ppval(pp, reshape (xi, szx));
     endif
   case "linear"
     dy = diff (y);
-    dx = diff (x);
-    if (pp)
-      coefs = [dy./dx, y(1:nx-1)];
-      xx = x;
-      if (have_jumps)
-        ## Omit zero-size intervals.
-        coefs(jumps) = [];
-        xx(jumps) = [];
-      endif
-      yi = mkpp (xx, coefs, szy(2:end));
+    dx = diff (x);    
+    dx = repmat (dx, [1 size(dy)(2:end)]);
+    coefs = [(dy./dx).'(:), y(1:nx-1, :).'(:)];
+    xx = x;
+
+    if (have_jumps)
+      ## Omit zero-size intervals.
+      coefs(jumps, :) = [];
+      xx(jumps) = [];
+    endif
+
+    pp = mkpp (xx, coefs, szy(2:end));
+    pp.orient = "first";
+
+    if (ispp)
+      yi = pp;
     else
-      ## find the interval containing the test point
-      idx = lookup (x, xi, "lr");
-      ## use the endpoints of the interval to define a line
-      s = (xi - x(idx))./dx(idx);
-      yi = bsxfun (@times, s, dy(idx,:)) + y(idx,:);
-      if (have_jumps)
-        ## Fix the corner cases of discontinuities at boundaries.
-        ## Internal discontinuities already handled correctly.
-        if (jumps (1))
-          mask = xi < x(1);
-          yi(mask,:) = y(1*ones (1, sum (mask)),:);
-        endif
-        if (jumps(nx-1))
-          mask = xi >= x(nx);
-          yi(mask,:) = y(nx*ones (1, sum (mask)),:);
-        endif
-      endif
+      yi = ppval(pp, reshape (xi, szx));
     endif
+
   case "*linear"
     dy = diff (y);
-    if (pp)
-      yi = mkpp (x(1) + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
+    coefs = [(dy/dx).'(:), y(1:nx-1, :).'(:)];
+    pp = mkpp (x, coefs, szy(2:end));
+    pp.orient = "first";
+
+    if (ispp)
+      yi = pp;
     else
-      ## find the interval containing the test point
-      t = (xi - x(1))/dx + 1;
-      idx = max (1, min (ny - 1, floor (t)));
+      yi = ppval(pp, reshape (xi, szx));
+    endif
 
-      ## use the endpoints of the interval to define a line
-      s = t - idx;
-      yi = bsxfun (@times, s, dy(idx,:)) + y(idx,:);
-    endif
-  case {"pchip", "*pchip"}
+  case {"pchip", "*pchip", "cubic", "*cubic"}
     if (nx == 2 || starmethod)
       x = linspace (x(1), x(nx), ny);
     endif
-    ## Note that pchip's arguments are transposed relative to interp1
-    if (pp)
-      yi = pchip (x.', y.');
-      yi.d = szy(2:end);
-    else
-      yi = pchip (x.', y.', xi.').';
-    endif
-
-  case {"cubic", "*cubic"}
-    if (nx < 4 || ny < 4)
-      error ("interp1: table too short");
-    endif
-
-    ## FIXME Is there a better way to treat pp return and *cubic
-    if (starmethod && ! pp)
-      ## From: Miloje Makivic
-      ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
-      t = (xi - x(1))/dx + 1;
-      idx = max (min (floor (t), ny-2), 2);
-      t = t - idx;
-      t2 = t.*t;
-      tp = 1 - 0.5*t;
-      a = (1 - t2).*tp;
-      b = (t2 + t).*tp;
-      c = (t2 - t).*tp/3;
-      d = (t2 - 1).*t/6;
-      J = ones (1, nc);
-
-      yi = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
-      + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
+    
+    if (ispp)
+      y = shiftdim (reshape (y, szy), 1);
+      yi = pchip (x, y);
     else
-      if (starmethod)
-        x = linspace (x(1), x(nx), ny).';
-        nx = ny;
-      endif
-
-      idx = lookup (x(2:nx-1), xi, "lr");
-
-      ## Construct cubic equations for each interval using divided
-      ## differences (computation of c and d don't use divided differences
-      ## but instead solve 2 equations for 2 unknowns). Perhaps
-      ## reformulating this as a lagrange polynomial would be more efficient.
-      i = 1:nx-3;
-      J = ones (1, nc);
-      dx = diff (x);
-      dx2 = x(i+1).^2 - x(i).^2;
-      dx3 = x(i+1).^3 - x(i).^3;
-      a = diff (y, 3)./dx(i,J).^3/6;
-      b = (diff (y(1:nx-1,:), 2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
-      c = (diff (y(1:nx-2,:), 1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
-      d = y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
-
-      if (pp)
-        xs = [x(1);x(3:nx-2)];
-        yi = mkpp ([x(1);x(3:nx-2);x(nx)],
-                   [a(:), (b(:) + 3.*xs(:,J).*a(:)), ...
-                    (c(:) + 2.*xs(:,J).*b(:) + 3.*xs(:,J)(:).^2.*a(:)), ...
-                    (d(:) + xs(:,J).*c(:) + xs(:,J).^2.*b(:) + ...
-                     xs(:,J).^3.*a(:))], szy(2:end));
-      else
-        yi = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
-              + c(idx,:)).*xi(:,J) + d(idx,:);
-      endif
+      y = shiftdim (y, 1);
+      yi = pchip (x, y, reshape (xi, szx));
     endif
   case {"spline", "*spline"}
     if (nx == 2 || starmethod)
       x = linspace(x(1), x(nx), ny);
     endif
-    ## Note that spline's arguments are transposed relative to interp1
-    if (pp)
-      yi = spline (x.', y.');
-      yi.d = szy(2:end);
+    
+    if (ispp)
+      y = shiftdim (reshape (y, szy), 1);
+      yi = spline (x, y);
     else
-      yi = spline (x.', y.', xi.').';
+      y = shiftdim (y, 1);
+      yi = spline (x, y, reshape (xi, szx));
     endif
   otherwise
     error ("interp1: invalid method '%s'", method);
   endswitch
 
-  if (! pp)
+  if (! ispp)
     if (! ischar (extrap))
       ## determine which values are out of range and set them to extrap,
       ## unless extrap == "extrap".
@@ -339,10 +276,24 @@
       maxx = max (x(1), x(nx));
 
       outliers = xi < minx | ! (xi <= maxx); # this catches even NaNs
-      yi(outliers, :) = extrap;
+      if (size_equal (outliers, yi))
+        yi(outliers) = extrap;
+        yi = reshape (yi, szx);
+      elseif (!isvector (yi))
+        if (strcmp (method, "pchip") || strcmp (method, "*pchip")
+          ||strcmp (method, "cubic") || strcmp (method, "*cubic")
+          ||strcmp (method, "spline") || strcmp (method, "*spline"))
+          yi(:, outliers) = extrap;
+          yi = shiftdim(yi, 1);
+        else
+          yi(outliers, :) = extrap;
+        endif
+      else
+        yi(outliers.') = extrap;
+      endif
     endif
-
-    yi = reshape (yi, [szx, szy(2:end)]);
+  else
+    yi.orient = "first";
   endif
 
 endfunction
@@ -394,6 +345,7 @@
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
+##FIXME: add test for n-d arguments here
 
 ## For each type of interpolated test, confirm that the interpolated
 ## value at the knots match the values at the knots.  Points away
@@ -595,7 +547,6 @@
 %!assert (interp1(1:2,1:2,1.4,"nearest"),1);
 %!error interp1(1,1,1, "linear");
 %!assert (interp1(1:2,1:2,1.4,"linear"),1.4);
-%!error interp1(1:3,1:3,1, "cubic");
 %!assert (interp1(1:4,1:4,1.4,"cubic"),1.4);
 %!assert (interp1(1:2,1:2,1.1, "spline"), 1.1);
 %!assert (interp1(1:3,1:3,1.4,"spline"),1.4);
@@ -604,7 +555,6 @@
 %!assert (interp1(1:2:4,1:2:4,1.4,"*nearest"),1);
 %!error interp1(1,1,1, "*linear");
 %!assert (interp1(1:2:4,1:2:4,[0,1,1.4,3,4],"*linear"),[NA,1,1.4,3,NA]);
-%!error interp1(1:3,1:3,1, "*cubic");
 %!assert (interp1(1:2:8,1:2:8,1.4,"*cubic"),1.4);
 %!assert (interp1(1:2,1:2,1.3, "*spline"), 1.3);
 %!assert (interp1(1:2:6,1:2:6,1.4,"*spline"),1.4);
@@ -612,5 +562,5 @@
 %!assert (interp1([3,2,1],[3,2,2],2.5),2.5)
 
 %!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap"), [-2,0.5,4,3,1.5])
-%!assert (interp1 ([4,4,3,2,0],[0,1,4,2,1],[1.5,4,4.5], "linear"), [0,1,NA])
+%!assert (interp1 ([4,4,3,2,0],[0,1,4,2,1],[1.5,4,4.5], "linear"), [1.75,1,NA])
 %!assert (interp1 (0:4, 2.5), 1.5)
--- a/scripts/general/isa.m
+++ b/scripts/general/isa.m
@@ -75,3 +75,22 @@
 %!assert (isa (uint16 (13), "numeric"), true)
 %!assert (isa (uint32 (13), "numeric"), true)
 %!assert (isa (uint64 (13), "numeric"), true)
+
+%!assert (isa (double (13), "double"));
+%!assert (isa (single (13), "single"));
+%!assert (isa (int8 (13), "int8"));
+%!assert (isa (int16 (13), "int16"));
+%!assert (isa (int32 (13), "int32"));
+%!assert (isa (int64 (13), "int64"));
+%!assert (isa (uint8 (13), "uint8"));
+%!assert (isa (uint16 (13), "uint16"));
+%!assert (isa (uint32 (13), "uint32"));
+%!assert (isa (uint64 (13), "uint64"));
+%!assert (isa ("string", "char"));
+%!assert (isa (true, "logical"));
+%!assert (isa (false, "logical"));
+%!assert (isa ({1, 2}, "cell"));
+%!test
+%! a.b = 1;
+%! assert (isa (a, "struct"));
+
--- a/scripts/general/iscolumn.m
+++ b/scripts/general/iscolumn.m
@@ -26,8 +26,6 @@
 
 function retval = iscolumn (x)
 
-  retval = false;
-
   if (nargin != 1)
     print_usage ();
   endif
--- a/scripts/general/isdir.m
+++ b/scripts/general/isdir.m
@@ -23,10 +23,11 @@
 ## @end deftypefn
 
 function retval = isdir (f)
-  if (nargin == 1)
-    ## Exist returns an integer but isdir should return a logical.
-    retval = (exist (f, "dir") == 7);
-  else
+  if (nargin != 1)
     print_usage ("isdir");
   endif
+
+  ## Exist returns an integer but isdir should return a logical.
+  retval = (exist (f, "dir") == 7);
+
 endfunction
--- a/scripts/general/isequal.m
+++ b/scripts/general/isequal.m
@@ -24,12 +24,12 @@
 
 function retval = isequal (x1, varargin)
 
-  if (nargin > 1)
-    retval = __isequal__ (false, x1, varargin{:});
-  else
+  if (nargin < 2)
     print_usage ();
   endif
 
+  retval = __isequal__ (false, x1, varargin{:});
+
 endfunction
 
 ## test size and shape
--- a/scripts/general/isequalwithequalnans.m
+++ b/scripts/general/isequalwithequalnans.m
@@ -25,12 +25,12 @@
 
 function retval = isequalwithequalnans (x1, varargin)
 
-  if (nargin > 1)
-    retval = __isequal__ (true, x1, varargin{:});
-  else
+  if (nargin < 2)
     print_usage ();
   endif
 
+  retval = __isequal__ (true, x1, varargin{:});
+
 endfunction
 
 ## test for equality
--- a/scripts/general/isrow.m
+++ b/scripts/general/isrow.m
@@ -26,8 +26,6 @@
 
 function retval = isrow (x)
 
-  retval = false;
-
   if (nargin != 1)
     print_usage ();
   endif
--- a/scripts/general/isscalar.m
+++ b/scripts/general/isscalar.m
@@ -26,12 +26,12 @@
 
 function retval = isscalar (x)
 
-  if (nargin == 1)
-    retval = numel (x) == 1;
-  else
+  if (nargin != 1)
     print_usage ();
   endif
 
+  retval = numel (x) == 1;
+
 endfunction
 
 %!assert(isscalar (1));
--- a/scripts/general/issquare.m
+++ b/scripts/general/issquare.m
@@ -28,42 +28,35 @@
 
 function retval = issquare (x)
 
-  if (nargin == 1)
-    if (ndims (x) == 2)
-      [r, c] = size (x);
-      retval = r == c;
-    else
-      retval = false;
-    endif
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  if (ndims (x) == 2)
+    [r, c] = size (x);
+    retval = r == c;
   else
-    print_usage ();
+    retval = false;
   endif
 
 endfunction
 
+%!assert(issquare ([]));
 %!assert(issquare (1));
-
 %!assert(!(issquare ([1, 2])));
-
-%!assert(issquare ([]));
-
 %!assert(issquare ([1, 2; 3, 4]));
-
-%!test
-%! assert(issquare ("t"));
-
+%!assert(!(issquare ([1, 2; 3, 4; 5, 6])));
+%!assert(!(issquare (ones (3,3,3))));
+%!assert(issquare ("t"));
 %!assert(!(issquare ("test")));
-
-%!test
-%! assert(issquare (["test"; "ing"; "1"; "2"]));
-
+%!assert(issquare (["test"; "ing"; "1"; "2"]));
 %!test
 %! s.a = 1;
 %! assert(issquare (s));
-
-%!assert(!(issquare ([1, 2; 3, 4; 5, 6])));
+%!assert(issquare ({1, 2; 3, 4}));
+%!assert(sparse (([1, 2; 3, 4])));
 
+%% Test input validation
 %!error issquare ();
-
 %!error issquare ([1, 2; 3, 4], 2);
 
--- a/scripts/general/isvector.m
+++ b/scripts/general/isvector.m
@@ -28,15 +28,13 @@
 
 function retval = isvector (x)
 
-  retval = 0;
-
-  if (nargin == 1)
-    sz = size (x);
-    retval = (ndims (x) == 2 && (sz(1) == 1 || sz(2) == 1));
-  else
+  if (nargin != 1)
     print_usage ();
   endif
 
+  sz = size (x);
+  retval = (ndims (x) == 2 && (sz(1) == 1 || sz(2) == 1));
+
 endfunction
 
 %!assert(isvector (1));
--- a/scripts/general/triplequad.m
+++ b/scripts/general/triplequad.m
@@ -37,22 +37,25 @@
 ##
 ## The optional argument @var{quadf} specifies which underlying integrator
 ## function to use.  Any choice but @code{quad} is available and the default
-## is @code{quadgk}.
+## is @code{quadcc}.
 ##
 ## Additional arguments, are passed directly to @var{f}.  To use the default
-## value for @var{tol} or @var{quadf} one may pass an empty matrix ([]).
+## value for @var{tol} or @var{quadf} one may pass ':' or an empty matrix ([]).
 ## @seealso{dblquad, quad, quadv, quadl, quadgk, quadcc, trapz}
 ## @end deftypefn
 
-function q = triplequad(f, xa, xb, ya, yb, za, zb, tol, quadf, varargin)
+function q = triplequad (f, xa, xb, ya, yb, za, zb, tol = 1e-6, quadf = @quadcc, varargin)
+
   if (nargin < 7)
     print_usage ();
   endif
-  if (nargin < 8 || isempty (tol))
+
+  ## Allow use of empty matrix ([]) to indicate default 
+  if (isempty (tol))
     tol = 1e-6;
   endif
-  if (nargin < 9 || isempty (quadf))
-    quadf = @quadgk;
+  if (isempty (quadf))
+    quadf = @quadcc;
   endif
 
   inner = @__triplequad_inner__;
@@ -61,7 +64,8 @@
     varargin = {};
   endif
 
-  q = dblquad(@(y, z) inner (y, z, f, xa, xb, tol, quadf, varargin{:}),ya, yb, za, zb, tol);
+  q = dblquad (@(y, z) inner (y, z, f, xa, xb, tol, quadf, varargin{:}), ya, yb, za, zb, tol);
+
 endfunction
 
 function q = __triplequad_inner__ (y, z, f, xa, xb, tol, quadf, varargin)
@@ -71,8 +75,11 @@
   endfor
 endfunction
 
-%% These tests are too expensive to run normally. Disable them
-% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadgk), pi ^ (3/2) * erf(1).^3, 1e-6)
-% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadl), pi ^ (3/2) * erf(1).^3, 1e-6)
-% !#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadv), pi ^ (3/2) * erf(1).^3, 1e-6)
+ 
+%!assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadcc), pi ^ (3/2) * erf(1).^3, 1e-6)
 
+%% These tests are too expensive to run normally (~30 sec each).  Disable them
+#%!assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadgk), pi ^ (3/2) * erf(1).^3, 1e-6)
+#%!#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadl), pi ^ (3/2) * erf(1).^3, 1e-6)
+#%!#assert (triplequad (@(x,y,z) exp(-x.^2 - y.^2 - z.^2) , -1, 1, -1, 1, -1, 1, [],  @quadv), pi ^ (3/2) * erf(1).^3, 1e-6)
+
--- a/scripts/help/__makeinfo__.m
+++ b/scripts/help/__makeinfo__.m
@@ -17,8 +17,8 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {[@var{retval}, @var{status}] =} __makeinfo__ (@var{text}, @var{output_type})
-## @deftypefnx {Function File} {[@var{retval}, @var{status}] =} __makeinfo__ (@var{text}, @var{output_type}, @var{see_also})
+## @deftypefn  {Function File} {[@var{retval}, @var{status}] =} __makeinfo__ (@var{text})
+## @deftypefnx {Function File} {[@var{retval}, @var{status}] =} __makeinfo__ (@var{text}, @var{output_type})
 ## Undocumented internal function.
 ## @end deftypefn
 
@@ -33,13 +33,6 @@
 ## @t{"plain text"}. If @var{output_type} is @t{"texinfo"}, the @t{@@seealso}
 ## macro is expanded, but otherwise the text is unaltered.
 ##
-## If the optional argument @var{see_also} is present, it is used to expand the
-## Octave specific @t{@@seealso} macro. This argument must be a function handle,
-## that accepts a cell array of strings as input argument (each elements of the
-## array corresponds to the arguments to the @t{@@seealso} macro), and return
-## the expanded string. If this argument is not given, the @t{@@seealso} macro
-## will be expanded to the text
-##
 ## @example
 ## See also: arg1, arg2@, ...
 ## @end example
@@ -60,7 +53,7 @@
 function [retval, status] = __makeinfo__ (text, output_type = "plain text", see_also = [])
 
   ## Check input
-  if (nargin == 0)
+  if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
 
@@ -72,72 +65,18 @@
     error ("__makeinfo__: second input argument must be a string");
   endif
 
-  ## Define the function which expands @seealso macro
-  if (isempty (see_also))
-    if (strcmpi (output_type, "plain text"))
-      see_also = @simple_see_also;
-    else
-      see_also = @simple_see_also_with_refs;
-    endif
-  endif
-
-  if (!isa (see_also, "function_handle"))
-    error ("__makeinfo__: third input argument must be the empty matrix, or a function handle");
-  endif
-
   ## It seems like makeinfo sometimes gets angry if the first character
   ## on a line is a space, so we remove these.
   text = strrep (text, "\n ", "\n");
 
   ## Handle @seealso macro
-  SEE_ALSO = "@seealso";
-  starts = strfind (text, SEE_ALSO);
-  for start = fliplr (starts)
-    if (start == 1 || (text (start-1) != "@"))
-      bracket_start = find (text (start:end) == "{", 1);
-      stop = find (text (start:end) == "}", 1);
-      if (!isempty (stop) && !isempty (bracket_start))
-        stop += start - 1;
-        bracket_start += start - 1;
-      else
-        bracket_start = start + length (SEE_ALSO);
-        stop = find (text (start:end) == "\n", 1);
-        if (isempty (stop))
-          stop = length (text);
-        else
-          stop += start - 1;
-        endif
-      endif
-      see_also_args = text (bracket_start+1:(stop-1));
-      see_also_args = strtrim (strsplit (see_also_args, ","));
-      expanded = see_also (see_also_args);
-      text = strcat (text (1:start-1), expanded, text (stop+1:end));
-    endif
-  endfor
-
+  if (strcmpi (output_type, "plain text")) 
+    text = regexprep (text, '@seealso *\{([^}]*)\}', "\nSee also: $1.\n\n");
+  else
+    text = regexprep (text, '@seealso *\{([^}]*)\}', "\nSee also: @ref{$1}.\n\n");
+  endif
   ## Handle @nospell macro
-  NOSPELL = "@nospell";
-  starts = strfind (text, NOSPELL);
-  for start = fliplr (starts)
-    if (start == 1 || (text (start-1) != "@"))
-      bracket_start = find (text (start:end) == "{", 1);
-      stop = find (text (start:end) == "}", 1);
-      if (!isempty (stop) && !isempty (bracket_start))
-        stop += start - 1;
-        bracket_start += start - 1;
-      else
-        bracket_start = start + length (NOSPELL);
-        stop = find (text (start:end) == "\n", 1);
-        if (isempty (stop))
-          stop = length (text);
-        else
-          stop += start - 1;
-        endif
-      endif
-      text(stop) = [];
-      text(start:bracket_start) = [];
-    endif
-  endfor
+  text = regexprep (text, '@nospell *\{([^}]*)\}', "$1");
 
   if (strcmpi (output_type, "texinfo"))
     status = 0;
@@ -180,12 +119,3 @@
   end_unwind_protect
 endfunction
 
-function expanded = simple_see_also (args)
-  expanded = strcat ("\nSee also:", sprintf (" %s,", args {:}));
-  expanded = strcat (expanded (1:end-1), "\n\n");
-endfunction
-
-function expanded = simple_see_also_with_refs (args)
-  expanded = strcat ("\nSee also:", sprintf (" @ref{%s},", args {:}));
-  expanded = strcat (expanded (1:end-1), "\n\n");
-endfunction
--- a/scripts/image/image.m
+++ b/scripts/image/image.m
@@ -19,9 +19,9 @@
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {} image (@var{img})
 ## @deftypefnx {Function File} {} image (@var{x}, @var{y}, @var{img})
-## Display a matrix as a color image.  The elements of @var{x} are indices
+## Display a matrix as a color image.  The elements of @var{img} are indices
 ## into the current colormap, and the colormap will be scaled so that the
-## extremes of @var{x} are mapped to the extremes of the colormap.
+## extremes of @var{img} are mapped to the extremes of the colormap.
 ##
 ## The axis values corresponding to the matrix elements are specified in
 ## @var{x} and @var{y}.  If you're not using gnuplot 4.2 or later, these
--- a/scripts/polynomial/mkpp.m
+++ b/scripts/polynomial/mkpp.m
@@ -17,50 +17,66 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {@var{pp} =} mkpp (@var{x}, @var{p})
-## @deftypefnx {Function File} {@var{pp} =} mkpp (@var{x}, @var{p}, @var{d})
+## @deftypefn  {Function File} {@var{pp} =} mkpp (@var{breaks}, @var{coefs})
+## @deftypefnx {Function File} {@var{pp} =} mkpp (@var{breaks}, @var{coefs}, @var{d})
+##
+## Construct a piece-wise polynomial (pp) structure from sample points
+## @var{breaks} and coefficients @var{coefs}.  @var{breaks} must be a vector of
+## strictly increasing values. The number of intervals is given by 
+## @code{@var{ni} = length (@var{breaks}) - 1}.
+## When @var{m} is the polynomial order @var{coefs} must be of 
+## size: @var{ni} x @var{m} + 1.
 ##
-## Construct a piecewise polynomial structure from sample points
-## @var{x} and coefficients @var{p}.  The i-th row of @var{p},
-## @code{@var{p} (@var{i},:)}, contains the coefficients for the polynomial
-## over the @var{i}-th interval, ordered from highest to
-## lowest.  There must be one row for each interval in @var{x}, so
-## @code{rows (@var{p}) == length (@var{x}) - 1}.
+## The i-th row of @var{coefs},
+## @code{@var{coefs} (@var{i},:)}, contains the coefficients for the polynomial
+## over the @var{i}-th interval, ordered from highest (@var{m}) to 
+## lowest (@var{0}).
 ##
-## @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}.
+## @var{coefs} may also be a multi-dimensional array, specifying a vector-valued
+## or array-valued polynomial. In that case the polynomial order is defined
+## by the length of the last dimension of @var{coefs}.
+## The size of first dimension(s) are given by the scalar or
+## vector @var{d}. If @var{d} is not given it is set to @code{1}. 
+## In any case @var{coefs} is reshaped to a 2d matrix of
+## size @code{[@var{ni}*prod(@var{d} @var{m})] }
 ##
 ## @seealso{unmkpp, ppval, spline}
 ## @end deftypefn
 
 function pp = mkpp (x, P, d)
+
+  # check number of arguments
   if (nargin < 2 || nargin > 3)
     print_usage ();
   endif
-  pp.x = x(:);
-  n = length (x) - 1;
-  if (n < 1)
+
+  # check x
+  if (length (x) < 2)
     error ("mkpp: at least one interval is needed");
   endif
-  nd = ndims (P);
-  k = size (P, nd);
-  if (nargin < 3)
-    if (nd == 2)
-      d = 1;
-    else
-      d = prod (size (P)(1:nd-1));
-    endif
+
+  if (!isvector (x))
+    error ("mkpp: x must be a vector");
   endif
-  pp.d = d;
-  pp.P = P = reshape (P, prod (d), [], k);
-  pp.orient = 0;
+
+  len = length (x) - 1;
+  dP = length (size (P));
 
-  if (size (P, 2) != n)
-    error ("mkpp: num intervals in X doesn't match num polynomials in P");
-  endif
+  pp = struct ("form", "pp",
+               "breaks", x(:).',
+               "coefs", [],
+               "pieces", len,
+               "order", prod (size (P)) / len,
+               "dim", 1);
+
+  if (nargin == 3)
+    pp.dim = d;
+    pp.order /= prod (d);
+  endif 
+
+  dim_vec = [pp.pieces * prod(pp.dim), pp.order];
+  pp.coefs = reshape (P, dim_vec);
+
 endfunction
 
 %!demo # linear interpolation
@@ -72,3 +88,25 @@
 %! xi=linspace(0,pi,50);
 %! plot(x,t,"x",xi,ppval(pp,xi));
 %! legend("control","interp");
+
+%!shared b,c,pp
+%! b = 1:3; c = 1:24; pp=mkpp(b,c);
+%!assert (pp.pieces,2);
+%!assert (pp.order,12);
+%!assert (pp.dim,1);
+%!assert (size(pp.coefs),[2,12]);
+%! pp=mkpp(b,c,2);
+%!assert (pp.pieces,2);
+%!assert (pp.order,6);
+%!assert (pp.dim,2);
+%!assert (size(pp.coefs),[4,6]);
+%! pp=mkpp(b,c,3);
+%!assert (pp.pieces,2);
+%!assert (pp.order,4);
+%!assert (pp.dim,3);
+%!assert (size(pp.coefs),[6,4]);
+%! pp=mkpp(b,c,[2,3]);
+%!assert (pp.pieces,2);
+%!assert (pp.order,2);
+%!assert (pp.dim,[2,3]);
+%!assert (size(pp.coefs),[12,2]);
--- a/scripts/polynomial/pchip.m
+++ b/scripts/polynomial/pchip.m
@@ -27,8 +27,8 @@
 ##
 ## The variable @var{x} must be a strictly monotonic vector (either
 ## increasing or decreasing).  While @var{y} can be either a vector or
-## array.  In the case where @var{y} is a vector, it must have a length
-## of @var{n}.  If @var{y} is an array, then the size of @var{y} must
+## an array.  In the case where @var{y} is a vector, it must have the
+## length @var{n}.  If @var{y} is an array, then the size of @var{y} must
 ## have the form
 ## @tex
 ## $$[s_1, s_2, \cdots, s_k, n]$$
@@ -73,15 +73,22 @@
     print_usage ();
   endif
 
+  ## make row vector
   x = x(:).';
   n = length (x);
 
   ## Check the size and shape of y
   if (isvector (y))
-    y = y(:).';
+    y = y(:).'; ##row vector
     szy = size (y);
+    if !(size_equal (x, y))
+      error ("pchip: length of X and Y must match")
+    endif
   else
     szy = size (y);
+    if (n != szy(end))
+      error ("pchip: length of X and last dimension of Y must match")
+    endif
     y = reshape (y, [prod(szy(1:end-1)), szy(end)]);
   endif
 
@@ -94,16 +101,12 @@
     error("pchip: X must be strictly monotonic");
   endif
 
-  if (columns (y) != n)
-    error("pchip: size of X and Y must match");
-  endif
-
-  f1 = y(:,1:n-1);
+  f1 = y(:, 1:n-1);
 
   ## Compute derivatives.
   d = __pchip_deriv__ (x, y, 2);
-  d1 = d(:,1:n-1);
-  d2 = d(:,2:n);
+  d1 = d(:, 1:n-1);
+  d2 = d(:, 2:n);
 
   ## This is taken from SLATEC.
   h = diag (h);
@@ -114,14 +117,12 @@
   c3 = del1 + del2;
   c2 = -c3 - del1;
   c3 = c3 / h;
-
   coeffs = cat (3, c3, c2, d1, f1);
-  pp = mkpp (x, coeffs, szy(1:end-1));
 
-  if (nargin == 2)
-    ret = pp;
-  else
-    ret = ppval (pp, xi);
+  ret = mkpp (x, coeffs, szy(1:end-1));
+
+  if (nargin == 3)
+    ret = ppval (ret, xi);
   endif
 
 endfunction
@@ -138,7 +139,7 @@
 %! %-------------------------------------------------------------------
 %! % confirm that pchip agreed better to discontinuous data than spline
 
-%!shared x,y
+%!shared x,y,y2,pp,yi1,yi2,yi3
 %! x = 0:8;
 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
 %!assert (pchip(x,y,x), y);
@@ -148,3 +149,23 @@
 %!assert (isempty(pchip(x',y',[])));
 %!assert (isempty(pchip(x,y,[])));
 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
+%!assert (pchip(x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
+%!assert (pchip(x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
+%!assert (pchip(x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
+%!test
+%! x=(0:8)*pi/4;y=[sin(x);cos(x)];
+%! y2(:,:,1)=y;y2(:,:,2)=y+1;y2(:,:,3)=y-1;
+%! pp=pchip(x,shiftdim(y2,2));
+%! yi1=ppval(pp,(1:4)*pi/4);
+%! yi2=ppval(pp,repmat((1:4)*pi/4,[5,1]));
+%! yi3=ppval(pp,[pi/2,pi]);
+%!assert(size(pp.coefs),[48,4]);
+%!assert(pp.pieces,8);
+%!assert(pp.order,4);
+%!assert(pp.dim,[3,2]);
+%!assert(ppval(pp,pi),[0,-1;1,0;-1,-2],1e-14);
+%!assert(yi3(:,:,2),ppval(pp,pi),1e-14);
+%!assert(yi3(:,:,1),[1,0;2,1;0,-1],1e-14);
+%!assert(squeeze(yi1(1,2,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
+%!assert(size(yi2),[3,2,5,4]);
+%!assert(squeeze(yi2(1,2,3,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
\ No newline at end of file
--- a/scripts/polynomial/ppder.m
+++ b/scripts/polynomial/ppder.m
@@ -17,28 +17,54 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{ppd} =} ppder (@var{pp})
-## Compute the piecewise derivative of the piecewise polynomial struct @var{pp}.
+## @deftypefn {Function File} {ppd =} ppder (pp, m)
+## Computes the piecewise @var{m}-th derivative of a piecewise polynomial
+## struct @var{pp}. If @var{m} is omitted the first derivate is
+## calculated.
 ## @seealso{mkpp, ppval, ppint}
 ## @end deftypefn
 
-function ppd = ppder (pp)
-  if (nargin != 1)
+function ppd = ppder (pp, m)
+
+  if ((nargin < 1) || (nargin > 2))
     print_usage ();
+  elseif (nargin == 1)
+    m = 1;
   endif
-  if (! isstruct (pp))
+
+  if !(isstruct (pp) && strcmp (pp.form, "pp"))
     error ("ppder: PP must be a structure");
   endif
 
   [x, p, n, k, d] = unmkpp (pp);
-  p = reshape (p, [], k);
-  if (k <= 1)
-    pd = zeros (rows (p), 1);
-    k = 1;
+
+  if (k - m <= 0)
+    x = [x(1) x(end)];
+    pd = zeros (prod (d), 1);
   else
-    k -= 1;
-    pd = p(:,1:k) * diag (k:-1:1);
+    f = k : -1 : 1;
+    ff = bincoeff (f, m + 1) .* factorial (m + 1) ./ f;
+    k -= m;
+    pd = p(:,1:k) * diag (ff(1:k));
   endif
+
   ppd = mkpp (x, pd, d);
 endfunction
 
+%!shared x,y,pp,ppd
+%! x=0:8;y=[x.^2;x.^3+1];pp=spline(x,y);
+%! ppd=ppder(pp);
+%!assert(ppval(ppd,x),[2*x;3*x.^2],1e-14)
+%!assert(ppd.order,3)
+%! ppd=ppder(pp,2);
+%!assert(ppval(ppd,x),[2*ones(size(x));6*x],1e-14)
+%!assert(ppd.order,2)
+%! ppd=ppder(pp,3);
+%!assert(ppd.order,1)
+%!assert(ppd.pieces,8)
+%!assert(size(ppd.coefs),[16,1])
+%! ppd=ppder(pp,4);
+%!assert(ppd.order,1)
+%!assert(ppd.pieces,1)
+%!assert(size(ppd.coefs),[2,1])
+%!assert(ppval(ppd,x),zeros(size(y)),1e-14)
--- a/scripts/polynomial/ppint.m
+++ b/scripts/polynomial/ppint.m
@@ -28,7 +28,7 @@
   if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
-  if (! isstruct (pp))
+  if (! isstruct (pp) && strcmp (pp.form, "pp"))
     error ("ppint: PP must be a structure");
   endif
 
@@ -39,17 +39,20 @@
   pi = p / diag (k:-1:1);
   k += 1;
   if (nargin == 1)
-    pi(:,k) = 0;
+    pi(:, k) = 0;
   else
-    pi(:,k) = repmat (c(:), n, 1);
+    pi(:, k) = repmat (c(:), n, 1);
   endif
 
   ppi = mkpp (x, pi, d);
 
-  ## Adjust constants so the the result is continuous.
-
-  jumps = reshape (ppjumps (ppi), prod (d), n-1);
-  ppi.P(:,2:n,k) -= cumsum (jumps, 2);
-
+  tmp = -cumsum (ppjumps (ppi), length (d) + 1);
+  ppi.coefs(prod(d)+1:end, k) = tmp(:);
+  
 endfunction
 
+%!shared x,y,pp,ppi
+%! x=0:8;y=[ones(size(x));x+1];pp=spline(x,y);
+%! ppi=ppint(pp);
+%!assert(ppval(ppi,x),[x;0.5*x.^2+x],1e-14)
+%!assert(ppi.order,5)
--- a/scripts/polynomial/ppjumps.m
+++ b/scripts/polynomial/ppjumps.m
@@ -28,29 +28,31 @@
   if (nargin != 1)
     print_usage ();
   endif
-  if (! isstruct (pp))
+  
+  if (! isstruct (pp) && strcmp (pp.form, "pp"))
     error ("ppjumps: PP must be a structure");
   endif
 
   ## Extract info.
-  x = pp.x;
-  P = pp.P;
-  d = pp.d;
-  [nd, n, k] = size (P);
+  [x, P, n, k, d] = unmkpp(pp);
+  nd = length (d) + 1;
 
   ## Offsets.
-  dx = diff (x(1:n)).';
-  dx = dx(ones (1, nd), :); # spread (do nothing in 1D)
+  dx = diff(x(1:n));
+  dx = repmat (dx, [prod(d), 1]);
+  dx = reshape (dx, [d, n-1]);
+  dx = shiftdim (dx, nd - 1);
 
-  ## Use Horner scheme to get limits from the left.
-  llim = P(:,1:n-1,1);
-  for i = 2:k;
+  ## Use Horner scheme.
+  if (k>1)
+    llim = shiftdim (reshape (P(1:(n-1) * prod(d), 1), [d, n-1]), nd - 1);
+  endif
+
+  for i = 2 : k;
     llim .*= dx;
-    llim += P(:,1:n-1,i);
+    llim += shiftdim (reshape (P(1:(n-1) * prod (d), i), [d, n-1]), nd - 1);
   endfor
-
-  rlim = P(:,2:n,k); # limits from the right
-  jumps = reshape (rlim - llim, [d, n-1]);
-
+  
+  rlim = shiftdim (ppval (pp, x(2:end-1)), nd - 1);
+  jumps = shiftdim (rlim - llim, 1);
 endfunction
-
--- a/scripts/polynomial/ppval.m
+++ b/scripts/polynomial/ppval.m
@@ -18,16 +18,18 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi})
-## Evaluate piecewise polynomial @var{pp} at the points @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
+## Evaluate piece-wise polynomial structure @var{pp} at the points @var{xi}.  
+## If @var{pp} describes a scalar polynomial function, the result is an
+## array of the same shape as @var{xi}.
+## Otherwise, the size of the result is @code{[pp.dim, length(@var{xi})]} if
+## @var{xi} is a vector, or @code{[pp.dim, size(@var{xi})]} if it is a
+## multi-dimensional array.  
+##
+##, 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}
+## @seealso{mkpp, unmkpp, spline, pchip, interp1}
 ## @end deftypefn
 
 function yi = ppval (pp, xi)
@@ -35,48 +37,85 @@
   if (nargin != 2)
     print_usage ();
   endif
-  if (! isstruct (pp))
-    error ("ppval: PP must be a structure");
+  if (! isstruct (pp) && strcmp (pp.form, "pp"))
+    error ("ppval: expects a pp-form structure");
   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)];
+  [x, P, n, k, d] = unmkpp (pp);
+  
+  ## dimension checks
+  sxi = size (xi);
+  if (isvector (xi))
+    xi = xi(:).';
   endif
+  
+  nd = length (d);
 
   ## Determine intervals.
-  xi = xi(:);
   xn = numel (xi);
-
   idx = lookup (x, xi, "lr");
 
+  P = reshape (P, [d, n * k]);
+  P = shiftdim (P, nd);
+  P = reshape (P, [n, k, d]);
+  Pidx = P(idx(:), :);#2d matrix size x: coefs*prod(d) y: prod(sxi)
+  
+  if (isvector(xi))
+    Pidx = reshape (Pidx, [xn, k, d]);
+    Pidx = shiftdim (Pidx, 1);
+    dimvec = [d, xn];
+  else
+    Pidx = reshape (Pidx, [sxi, k, d]);
+    Pidx = shiftdim (Pidx, length (sxi));
+    dimvec = [d, sxi];
+  end
+  ndv = length (dimvec);
+
   ## Offsets.
-  dx = (xi - x(idx)).';
-  dx = dx(ones (1, nd), :); # spread (do nothing in 1D)
+  dx = (xi - x(idx));
+  dx = repmat (dx, [prod(d), 1]);
+  dx = reshape (dx, dimvec);
+  dx = shiftdim (dx, ndv - 1);
 
   ## Use Horner scheme.
-  yi = P(:,idx,1);
-  for i = 2:k;
+  yi = Pidx;
+  if (k > 1)
+    yi = shiftdim (reshape (Pidx(1,:), dimvec), ndv - 1);
+  endif
+  
+  for i = 2 : k;
     yi .*= dx;
-    yi += P(:,idx,i);
+    yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1);
   endfor
-
+  
   ## Adjust shape.
-  yi = reshape (yi, yisz);
-  if (d != 1 && pp.orient == 1)
-    ## Switch dimensions to match interp1 order.
-    yi = shiftdim (yi, length (d));
+  if ((numel (xi) > 1) || (length (d) == 1))
+    yi = reshape (shiftdim (yi, 1), dimvec);
   endif
 
+  if (isvector (xi) && (d == 1))
+    yi = reshape (yi, sxi);
+  elseif (isfield (pp, "orient") && strcmp (pp.orient, "first"))
+    yi = shiftdim(yi, nd);
+  endif
+
+  ##
+  #if (d == 1)
+  #  yi = reshape (yi, sxi);
+  #endif
+  
 endfunction
+
+%!shared b,c,pp,pp2,xi,abserr
+%! b = 1:3; c = ones(2); pp=mkpp(b,c);abserr = 1e-14;pp2=mkpp(b,[c;c],2);
+%! xi = [1.1 1.3 1.9 2.1];
+%!assert (ppval(pp,1.1), 1.1, abserr);
+%!assert (ppval(pp,2.1), 1.1, abserr);
+%!assert (ppval(pp,xi), [1.1 1.3 1.9 1.1], abserr);
+%!assert (ppval(pp,xi.'), [1.1 1.3 1.9 1.1].', abserr);
+%!assert (ppval(pp2,1.1), [1.1;1.1], abserr);
+%!assert (ppval(pp2,2.1), [1.1;1.1], abserr);
+%!assert (ppval(pp2,xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
+%!assert (ppval(pp2,xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
+%!assert (size(ppval(pp2,[xi;xi])), [2 2 4]);
--- a/scripts/polynomial/spline.m
+++ b/scripts/polynomial/spline.m
@@ -83,15 +83,15 @@
   ## Check the size and shape of y
   ndy = ndims (y);
   szy = size (y);
-  if (ndy == 2 && (szy(1) == 1 || szy(2) == 1))
-    if (szy(1) == 1)
+  if (ndy == 2 && (szy(1) == n || szy(2) == n))
+    if (szy(2) == n)
       a = y.';
     else
       a = y;
       szy = fliplr (szy);
     endif
   else
-    a = reshape (y, [prod(szy(1:end-1)), szy(end)]).';
+    a = shiftdim (reshape (y, [prod(szy(1:end-1)), szy(end)]), 1);
   endif
 
   for k = (1:columns (a))(any (isnan (a)))
@@ -120,9 +120,9 @@
 
     if (n == 2)
       d = (dfs + dfe) / (x(2) - x(1)) ^ 2 + ...
-	2 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 3;
+          2 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 3;
       c = (-2 * dfs - dfe) / (x(2) - x(1)) - ...
-	3 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 2;
+          3 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 2;
       b = dfs;
       a = a(1,:);
 
@@ -132,7 +132,7 @@
       a = a(1:n-1,:);
     else
       if (n == 3)
-	dg = 1.5 * h(1) - 0.5 * h(2);
+        dg = 1.5 * h(1) - 0.5 * h(2);
         c(2:n-1,:) = 1/dg(1);
       else
         dg = 2 * (h(1:n-2) .+ h(2:n-1));
@@ -153,9 +153,9 @@
       endif
 
       c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
-		- c(2,:) * h(1)) / (2 * h(1));
+             - c(2,:) * h(1)) / (2 * h(1));
       c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
-		  + c(n-1,:) * h(n-1)) / (2 * h(n-1));
+             + c(n-1,:) * h(n-1)) / (2 * h(n-1));
       b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
         - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
       d = diff (c) ./ (3 * h(1:n-1, idx));
@@ -229,15 +229,14 @@
           - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
       d = diff (c) ./ (3 * h(1:n-1, idx));
 
-      d = d(1:n-1,:);
-      c = c(1:n-1,:);
-      b = b(1:n-1,:);
-      a = a(1:n-1,:);
+      d = d(1:n-1,:);d = d.'(:);
+      c = c(1:n-1,:);c = c.'(:);
+      b = b(1:n-1,:);b = b.'(:);
+      a = a(1:n-1,:);a = a.'(:);
     endif
 
   endif
-  coeffs = cat (3, d.', c.', b.', a.');
-  ret = mkpp (x, coeffs, szy(1:end-1));
+  ret = mkpp (x, cat (2, d, c, b, a), szy(1:end-1));
 
   if (nargin == 3)
     ret = ppval (ret, xi);
@@ -263,6 +262,9 @@
 %!assert (isempty(spline(x',y',[])));
 %!assert (isempty(spline(x,y,[])));
 %!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
+%!assert (spline(x,[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
+%!assert (spline(x',[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
+%!assert (spline(x',[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
 %! y = cos(x) + i*sin(x);
 %!assert (spline(x,y,x), y, abserr)
 %!assert (real(spline(x,y,x)), real(y), abserr);
--- a/scripts/polynomial/unmkpp.m
+++ b/scripts/polynomial/unmkpp.m
@@ -50,15 +50,13 @@
   if (nargin == 0)
     print_usage ();
   endif
-  if (! isstruct (pp))
+  if (! (isstruct (pp) && strcmp (pp.form, "pp")))
     error ("unmkpp: expecting piecewise polynomial structure");
   endif
-  x = pp.x;
-  P = pp.P;
-  n = size (P, 2);
-  k = size (P, 3);
-  d = pp.d;
-  if (d == 1)
-    P = reshape (P, n, k);
-  endif
+  x = pp.breaks;
+  P = pp.coefs;
+  n = pp.pieces;
+  k = pp.order;
+  d = pp.dim;
+
 endfunction
--- a/scripts/statistics/base/center.m
+++ b/scripts/statistics/base/center.m
@@ -23,7 +23,7 @@
 ## If @var{x} is a vector, subtract its mean.
 ## If @var{x} is a matrix, do the above for each column.
 ## If the optional argument @var{dim} is given, operate along this dimension.
-## @seealso{studentize}
+## @seealso{zscore}
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
--- a/scripts/statistics/base/mode.m
+++ b/scripts/statistics/base/mode.m
@@ -90,7 +90,7 @@
   if (issparse (x))
     m = sparse (sz2(1), sz2(2));
   else
-    m = zeros (sz2);
+    m = zeros (sz2, class (x));
   endif
   for i = 1 : prod (sz2)
     c{i} = xs (t2 (:, i) == f(i), i);
--- a/scripts/statistics/base/module.mk
+++ b/scripts/statistics/base/module.mk
@@ -33,9 +33,9 @@
   statistics/base/spearman.m \
   statistics/base/statistics.m \
   statistics/base/std.m \
-  statistics/base/studentize.m \
   statistics/base/table.m \
-  statistics/base/var.m
+  statistics/base/var.m \
+  statistics/base/zscore.m
 
 FCN_FILES += $(statistics_base_FCN_FILES)
 
--- a/scripts/statistics/base/quantile.m
+++ b/scripts/statistics/base/quantile.m
@@ -308,6 +308,10 @@
     error ("quantile: X must be a numeric vector or matrix");
   endif
 
+  if (isinteger (x))
+    x = double (x);
+  endif
+
   ## Save length and set shape of quantiles.
   n = numel (p);
   p = p(:);
@@ -320,7 +324,7 @@
   nx = size (x, 2);
 
   ## Initialize output values.
-  inv = Inf*(-(p < 0) + (p > 1));
+  inv = Inf(class (x)) * (-(p < 0) + (p > 1));
   inv = repmat (inv, 1, nx);
 
   ## Do the work.
copy from scripts/statistics/base/studentize.m
copy to scripts/statistics/base/zscore.m
--- a/scripts/statistics/base/studentize.m
+++ b/scripts/statistics/base/zscore.m
@@ -17,8 +17,8 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {} studentize (@var{x})
-## @deftypefnx {Function File} {} studentize (@var{x}, @var{dim})
+## @deftypefn  {Function File} {} zscore (@var{x})
+## @deftypefnx {Function File} {} zscore (@var{x}, @var{dim})
 ## If @var{x} is a vector, subtract its mean and divide by its standard
 ## deviation.
 ##
@@ -31,14 +31,14 @@
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
 ## Description: Subtract mean and divide by standard deviation
 
-function t = studentize (x, dim)
+function t = zscore (x, dim)
 
   if (nargin != 1 && nargin != 2)
     print_usage ();
   endif
 
   if (! isnumeric(x))
-    error ("studentize: X must be a numeric vector or matrix");
+    error ("zscore: X must be a numeric vector or matrix");
   endif
 
   if (isinteger (x))
@@ -56,7 +56,7 @@
   else
     if (!(isscalar (dim) && dim == fix (dim))
         || !(1 <= dim && dim <= nd))
-      error ("studentize: DIM must be an integer and a valid dimension");
+      error ("zscore: DIM must be an integer and a valid dimension");
     endif
   endif
 
@@ -72,17 +72,17 @@
 
 endfunction
 
-%!assert(studentize ([1,2,3]), [-1,0,1])
-%!assert(studentize (int8 ([1,2,3])), [-1,0,1])
-#%!assert(studentize (ones (3,2,0,2)), zeros (3,2,0,2))
-%!assert(studentize ([2,0,-2;0,2,0;-2,-2,2]), [1,0,-1;0,1,0;-1,-1,1])
+%!assert(zscore ([1,2,3]), [-1,0,1])
+%!assert(zscore (int8 ([1,2,3])), [-1,0,1])
+#%!assert(zscore (ones (3,2,0,2)), zeros (3,2,0,2))
+%!assert(zscore ([2,0,-2;0,2,0;-2,-2,2]), [1,0,-1;0,1,0;-1,-1,1])
 
 %% Test input validation
-%!error studentize ()
-%!error studentize (1, 2, 3)
-%!error studentize ([true true])
-%!error studentize (1, ones(2,2))
-%!error studentize (1, 1.5)
-%!error studentize (1, 0)
-%!error studentize (1, 3)
+%!error zscore ()
+%!error zscore (1, 2, 3)
+%!error zscore ([true true])
+%!error zscore (1, ones(2,2))
+%!error zscore (1, 1.5)
+%!error zscore (1, 0)
+%!error zscore (1, 3)
 
--- a/scripts/strings/mat2str.m
+++ b/scripts/strings/mat2str.m
@@ -28,7 +28,7 @@
 ## scalar then both real and imaginary parts of the matrix are printed
 ## to the same precision.  Otherwise @code{@var{n} (1)} defines the
 ## precision of the real part and @code{@var{n} (2)} defines the
-## precision of the imaginary part.  The default for @var{n} is 17.
+## precision of the imaginary part.  The default for @var{n} is 15.
 ##
 ## If the argument 'class' is given, then the class of @var{x} is
 ## included in the string in such a way that the eval will result in the
@@ -56,13 +56,13 @@
 
   if (nargin < 2 || isempty (n))
     ## Default precision
-    n = 17;
+    n = 15;
   endif
 
   if (nargin < 3)
     if (ischar (n))
       cls = n;
-      n = 17;
+      n = 15;
     else
       cls = "";
     endif
@@ -137,3 +137,4 @@
 %!assert (mat2str (true), "true");
 %!assert (mat2str (false), "false");
 %!assert (mat2str (logical (eye (2))), "[true,false;false,true]");
+%!assert (mat2str (0.7), "0.7")
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,10 +1,15 @@
+2011-04-19  Kai Habel  <kai.habel@gmx.de>
+
+	* src/DLD-FUNCTIONS/__init_fltk__.cc(plot_window::plot_window):
+	Instantiate canvas before uimenu.
+
 2011-04-13  Rik  <octave@nomad.inbox5.com>
 
 	* help.cc: Add spaces after commas in @seealso blocks.
 
 2011-04-12  Rik  <octave@nomad.inbox5.com>
 
-	* load-path.cc (restoredefaultpath): Correct use of it's -> its in 
+	* load-path.cc (restoredefaultpath): Correct use of it's -> its in
 	documentation.
 
 2011-04-10  John Eaton  <jwe@octave.org>
@@ -90,7 +95,7 @@
 
 	* DLD-FUNCTIONS/inv.cc (inv, inverse), DLD-FUNCTIONS/tril.cc (tril),
 	data.cc (cumsum, szie), file-io.cc (fgets), ov-typeinfo.cc (typeinfo),
-	ov-usr-fcn.cc (nargout), utils.cc (make_absolute_filename), 
+	ov-usr-fcn.cc (nargout), utils.cc (make_absolute_filename),
 	variables.cc (who): Improve docstrings
 
 2011-03-25  John W. Eaton  <jwe@octave.org>
--- a/src/DLD-FUNCTIONS/__init_fltk__.cc
+++ b/src/DLD-FUNCTIONS/__init_fltk__.cc
@@ -648,13 +648,13 @@
     begin ();
     {
 
+      canvas = new
+        OpenGL_fltk (0, 0, ww , hh - status_h, number ());
+
       uimenu = new
         fltk_uimenu(0, 0, ww, menu_h);
       uimenu->hide ();
 
-      canvas = new
-        OpenGL_fltk (0, 0, ww , hh - status_h, number ());
-
       bottom = new
         Fl_Box (0,
                 hh - status_h,
--- a/src/DLD-FUNCTIONS/quadcc.cc
+++ b/src/DLD-FUNCTIONS/quadcc.cc
@@ -50,8 +50,7 @@
   int depth, rdepth, ndiv;
 } cquad_ival;
 
-/* Some constants and matrices that we'll need.
-    */
+/* Some constants and matrices that we'll need.  */
 
 static const double xi[33] = {
   -1., -0.99518472667219688624, -0.98078528040323044912,
@@ -1473,9 +1472,7 @@
 }
 
 
-
-/* The actual integration routine.
-    */
+/* The actual integration routine.  */
 
 DEFUN_DLD (quadcc, args, nargout,
 "-*- texinfo -*-\n\
@@ -1545,6 +1542,7 @@
 @seealso{quad, quadv, quadl, quadgk, trapz, dblquad, triplequad}\n\
 @end deftypefn")
 {
+  octave_value_list retval;
 
   /* Some constants that we will need. */
   static const int n[4] = { 4, 8, 16, 32 };
@@ -1563,11 +1561,11 @@
   double a, b, tol, iivals[cquad_heapsize], *sing;
 
   /* Variables needed for transforming the integrand. */
-  int wrap = 0;
+  bool wrap = false;
   double xw;
 
   /* Stuff we will need to call the integrand. */
-  octave_value_list fargs, retval;
+  octave_value_list fargs, fvals;
 
   /* Actual variables (as opposed to constants above). */
   double m, h, ml, hl, mr, hr, temp;
@@ -1580,48 +1578,49 @@
 
 
   /* Parse the input arguments. */
-  if (nargin < 1)
+  if (nargin < 3)
     {
-      error
-        ("quadcc: first argument (integrand) of type function handle required");
-      return octave_value_list ();
+      print_usage ();
+      return retval;
     }
+
+  if (args(0).is_function_handle () || args(0).is_inline_function ())
+    fcn = args(0).function_value ();
   else
     {
-      if (args (0).is_function_handle () || args (0).is_inline_function ())
-        fcn = args (0).function_value ();
-      else
-        {
-          error ("quadcc: first argument (integrand) must be a function handle or an inline function");
-          return octave_value_list();
-        }
+       std::string fcn_name = unique_symbol_name ("__quadcc_fcn_");
+       std::string fname = "function y = ";
+       fname.append (fcn_name);
+       fname.append ("(x) y = ");
+       fcn = extract_function (args(0), "quadcc", fcn_name, fname,
+                               "; endfunction");
     }
 
-  if (nargin < 2 || !args (1).is_real_scalar ())
+  if (!args(1).is_real_scalar ())
     {
-      error ("quadcc: second argument (left interval edge) must be a single real scalar");
-      return octave_value_list ();
+      error ("quadcc: lower limit of integration (A) must be a single real scalar");
+      return retval;
     }
   else
-    a = args (1).double_value ();
+    a = args(1).double_value ();
 
-  if (nargin < 3 || !args (2).is_real_scalar ())
+  if (!args(2).is_real_scalar ())
     {
-      error ("quadcc: third argument (right interval edge) must be a single real scalar");
-      return octave_value_list ();
+      error ("quadcc: upper limit of integration (B) must be a single real scalar");
+      return retval;
     }
   else
-    b = args (2).double_value ();
+    b = args(2).double_value ();
 
-  if (nargin < 4)
+  if (nargin < 4 || args(3).is_empty ())
     tol = 1.0e-6;
-  else if (!args (3).is_real_scalar ())
+  else if (!args(3).is_real_scalar () || args(3).double_value () <= 0)
     {
-      error ("quadcc: fourth argument (tolerance) must be a single real scalar");
-      return octave_value_list ();
+      error ("quadcc: tolerance (TOL) must be a single real scalar > 0");
+      return retval;
     }
   else
-    tol = args (3).double_value ();
+    tol = args(3).double_value ();
 
   if (nargin < 5)
     {
@@ -1629,20 +1628,21 @@
       iivals[0] = a;
       iivals[1] = b;
     }
-  else if (!(args (4).is_real_scalar () || args (4).is_real_matrix ()))
+  else if (!(args(4).is_real_scalar () || args(4).is_real_matrix ()))
     {
-      error ("quadcc: fifth argument (singularities) must be a vector of real values");
-      return octave_value_list ();
+      error ("quadcc: list of singularities (SING) must be a vector of real values");
+      return retval;
     }
   else
     {
-      nivals = 1 + args (4).length ();
-      if ( nivals > cquad_heapsize ) {
-        error ("quadcc: maximum number of singular points is limited to %i",
-               cquad_heapsize-1);
-        return octave_value_list();
+      nivals = 1 + args(4).length ();
+      if (nivals > cquad_heapsize) 
+        {
+          error ("quadcc: maximum number of singular points is limited to %i",
+                 cquad_heapsize-1);
+          return retval;
         }
-      sing = args (4).array_value ().fortran_vec ();
+      sing = args(4).array_value ().fortran_vec ();
       iivals[0] = a;
       for (i = 0; i < nivals - 2; i++)
         iivals[i + 1] = sing[i];
@@ -1652,7 +1652,7 @@
   /* If a or b are +/-Inf, transform the integral. */
   if (xisinf (a) || xisinf (b))
     {
-      wrap = 1;
+      wrap = true;
       for (i = 0; i <= nivals; i++)
         if (xisinf (iivals[i]))
           iivals[i] = copysign (1.0, iivals[i]);
@@ -1688,19 +1688,18 @@
           for (i = 0; i <= n[3]; i++)
             ex (i) = m + xi[i] * h;
         }
-      fargs (0) = ex;
-      retval = feval (fcn, fargs, 1);
-      if (retval.length () != 1 || !retval (0).is_real_matrix ())
+      fargs(0) = ex;
+      fvals = feval (fcn, fargs, 1);
+      if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
         {
-          error
-            ("quadcc: integrand must return a single, real-valued vector");
-          return octave_value_list ();
+          error ("quadcc: integrand F must return a single, real-valued vector");
+          return retval;
         }
-      Matrix effex = retval (0).matrix_value ();
+      Matrix effex = fvals(0).matrix_value ();
       if (effex.length () != ex.length ())
         {
-          error ("quadcc: integrand must return a single, real-valued vector of the same size as the input");
-          return octave_value_list ();
+          error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
+          return retval;
         }
       for (i = 0; i <= n[3]; i++)
         {
@@ -1809,18 +1808,18 @@
                 for (i = 0; i < n[d] / 2; i++)
                   ex (i) = m + xi[(2 * i + 1) * skip[d]] * h;
               }
-            fargs (0) = ex;
-            retval = feval (fcn, fargs, 1);
-            if (retval.length () != 1 || !retval (0).is_real_matrix ())
+            fargs(0) = ex;
+            fvals = feval (fcn, fargs, 1);
+            if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector");
+                return retval;
               }
-            Matrix effex = retval (0).matrix_value ();
+            Matrix effex = fvals(0).matrix_value ();
             if (effex.length () != ex.length ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector of the same size as the input");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
+                return retval;
               }
             neval += effex.length ();
             for (i = 0; i < n[d] / 2; i++)
@@ -1957,18 +1956,18 @@
                 for (i = 0; i < n[0] - 1; i++)
                   ex (i) = ml + xi[(i + 1) * skip[0]] * hl;
               }
-            fargs (0) = ex;
-            retval = feval (fcn, fargs, 1);
-            if (retval.length () != 1 || !retval (0).is_real_matrix ())
+            fargs(0) = ex;
+            fvals = feval (fcn, fargs, 1);
+            if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector");
+                return retval;
               }
-            Matrix effex = retval (0).matrix_value ();
+            Matrix effex = fvals(0).matrix_value ();
             if (effex.length () != ex.length ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector of the same size as the input");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
+                return retval;
               }
             neval += effex.length ();
             for (i = 0; i < n[0] - 1; i++)
@@ -2053,18 +2052,18 @@
                 for (i = 0; i < n[0] - 1; i++)
                   ex (i) = mr + xi[(i + 1) * skip[0]] * hr;
               }
-            fargs (0) = ex;
-            retval = feval (fcn, fargs, 1);
-            if (retval.length () != 1 || !retval (0).is_real_matrix ())
+            fargs(0) = ex;
+            fvals = feval (fcn, fargs, 1);
+            if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector");
+                return retval;
               }
-            Matrix effex = retval (0).matrix_value ();
+            Matrix effex = fvals(0).matrix_value ();
             if (effex.length () != ex.length ())
               {
-                error ("quadcc: integrand must return a single, real-valued vector of the same size as the input");
-                return octave_value_list ();
+                error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
+                return retval;
               }
             neval += effex.length ();
             for (i = 0; i < n[0] - 1; i++)
@@ -2234,11 +2233,39 @@
     }
 */
   /* Clean up and present the results. */
-  retval (0) = igral;
+  if (nargout > 2)
+    retval(2) = neval;
   if (nargout > 1)
-    retval (1) = err;
-  if (nargout > 2)
-    retval (2) = neval;
+    retval(1) = err;
+  retval(0) = igral;
   /* All is well that ends well. */
   return retval;
 }
+
+
+/* 
+
+%!assert (quadcc(@sin,-pi,pi), 0, 1e-6)
+%!assert (quadcc(inline('sin'),-pi,pi), 0, 1e-6)
+%!assert (quadcc('sin',-pi,pi), 0, 1e-6)
+
+%!assert (quadcc(@sin,-pi,0), -2, 1e-6)
+%!assert (quadcc(@sin,0,pi), 2, 1e-6)
+%!assert (quadcc(@(x) 1./sqrt(x), 0, 1), 2, 1e-6)
+%!assert (quadcc(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6)
+
+%!assert (quadcc (@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi), 1e-6)
+%!assert (quadcc (@(x) exp(-x .^ 2), -Inf, 0), sqrt(pi)/2, 1e-6)
+
+%% Test input validation
+%!error (quadcc ())
+%!error (quadcc (@sin))
+%!error (quadcc (@sin, 0))
+%!error (quadcc (@sin, ones(2), pi))
+%!error (quadcc (@sin, -i, pi))
+%!error (quadcc (@sin, 0, ones(2)))
+%!error (quadcc (@sin, 0, i))
+%!error (quadcc (@sin, 0, pi, 0))
+%!error (quadcc (@sin, 0, pi, 1e-6, [ i ]))
+
+*/
--- a/test/ChangeLog
+++ b/test/ChangeLog
@@ -1,3 +1,9 @@
+2011-04-11  Rik  <octave@nomad.inbox5.com>
+
+	* fntests.m: Remove deprecated and private functions from list of 
+	functions requiring tests.  Count functions with %!demo blocks as
+	having tests.
+
 2011-04-03  Rik  <octave@nomad.inbox5.com>
 
 	* test_diag_perm.m: Reverse previous changeset.  Return 3-input form
--- a/test/fntests.m
+++ b/test/fntests.m
@@ -97,6 +97,17 @@
   endif
 endfunction
 
+function retval = has_demos (f)
+  fid = fopen (f);
+  if (fid >= 0)
+    str = fread (fid, "*char")';
+    fclose (fid);
+    retval = ! isempty (regexp (str, '^%!demo', "lineanchors"));
+  else
+    error ("fopen failed: %s", f);
+  endif
+endfunction
+
 function [dp, dn, dxf, dsk] = run_test_dir (fid, d);
   global files_with_tests;
   global files_with_no_tests;
@@ -113,6 +124,8 @@
         [p, n, xf, sk] = test (nm(1:(end-2)), "quiet", fid);
         print_pass_fail (n, p);
         files_with_tests(end+1) = ffnm;
+      elseif (has_demos (ffnm))
+        files_with_tests(end+1) = ffnm;
       else
         files_with_no_tests(end+1) = ffnm;
       endif
@@ -164,6 +177,8 @@
         dxf += xf;
         dsk += sk;
         files_with_tests(end+1) = f;
+      elseif (has_demos (f))
+        files_with_tests(end+1) = f;
       elseif (has_functions (f))
         ## To reduce the list length, only mark .cc files that contain
         ## DEFUN definitions.
@@ -192,20 +207,15 @@
 endfunction
 
 function n = num_elts_matching_pattern (lst, pat)
-  n = 0;
-  for i = 1:length (lst)
-    if (! isempty (regexp (lst{i}, pat, "once")))
-      n++;
-    endif
-  endfor
+  n = sum (cellfun (@(x) !isempty (x), regexp (lst, pat, 'once')));
 endfunction
 
 function report_files_with_no_tests (with, without, typ)
-  pat = cstrcat ("\\", typ, "$");
+  pat = cstrcat ('\', typ, "$");
   n_with = num_elts_matching_pattern (with, pat);
   n_without = num_elts_matching_pattern (without, pat);
   n_tot = n_with + n_without;
-  printf ("\n%d (of %d) %s files have no tests.\n", n_without, n_tot, typ);
+  printf ("\n%d (of %d) %s files have no tests or demos.\n", n_without, n_tot, typ);
 endfunction
 
 pso = page_screen_output ();
@@ -258,6 +268,12 @@
     puts ("because the needed libraries were not present when Octave was built.\n");
   endif
 
+  ## Weed out deprecated and private functions
+  weed_idx = cellfun (@isempty, regexp (files_with_tests, '\bdeprecated\b|\bprivate\b', 'once'));
+  files_with_tests = files_with_tests(weed_idx);
+  weed_idx = cellfun (@isempty, regexp (files_with_no_tests, '\bdeprecated\b|\bprivate\b', 'once'));
+  files_with_no_tests = files_with_no_tests(weed_idx);
+
   report_files_with_no_tests (files_with_tests, files_with_no_tests, ".m");
   report_files_with_no_tests (files_with_tests, files_with_no_tests, ".cc");