--- a/libcruft/ChangeLog
+++ b/libcruft/ChangeLog
@@ -1,3 +1,7 @@
+2006-05-22  John W. Eaton  <jwe@octave.org>
+
+	* lapack/dlantr.f, lapack/zlantr.f: New files.
+
 2006-05-03  David Bateman  <dbateman@free.fr>
 
 	* lapack/dpocon.f, lapack/zpocon.f, lapack/dpotrs.f, 

--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,19 @@
+2006-05-22  John W. Eaton  <jwe@octave.org>
+
+	* scripts/general/lookup.m: New file from Octave Forge.
+
+2006-05-22  David Bateman  <dbateman@free.fr>
+
+	* scripts/general/nthroot.m, scripts/linear-algebra/rref.m,
+	scripts/specfun/isprime.m, scripts/specfun/primes.m,
+	scripts/specfun/factor.m, scripts/specfun/nchoosek.m,
+	scripts/specfun/legendre.m, scripts/specfun/perms.m,
+	scripts/special-matrix/rosser.m,
+	scripts/special-matrix/wilkinson.m,
+	scripts/special-matrix/pascal.m, scripts/special-matrix/magic.m,
+	scripts/special-matrix/hadamard.m, scripts/strings/strtok.m:
+	New files from Octave Forge.
+
 2006-05-19  David Bateman  <dbateman@free.fr>
 
 	* polynomial/unmkpp.m, polynomial/mkpp.m, polynomial/spline.m,

new file mode 100644
--- /dev/null
+++ b/scripts/general/lookup.m
@@ -0,0 +1,101 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{idx} =} lookup (@var{table}, @var{y})
+## Lookup values in a sorted table.  Usually used as a prelude to
+## interpolation.
+##
+## If table is strictly increasing and @code{idx = lookup (table, y)}, then
+## @code{table(idx(i)) <= y(i) < table(idx(i+1))} for all @code{y(i)}
+## within the table.  If @code{y(i)} is before the table, then 
+## @code{idx(i)} is 0. If @code{y(i)} is after the table then
+## @code{idx(i)} is @code{table(n)}.
+##
+## If the table is strictly decreasing, then the tests are reversed.
+## There are no guarantees for tables which are non-monotonic or are not
+## strictly monotonic.
+##
+## To get an index value which lies within an interval of the table,
+## use:
+##
+## @example
+## idx = lookup (table(2:length(table)-1), y) + 1
+## @end example
+##
+## @noindent
+## This expression puts values before the table into the first
+## interval, and values after the table into the last interval.
+## @end deftypefn
+
+## Changed from binary search to sort.
+## Thanks to Kai Habel <kai.habel@gmx.de> for the suggestion.
+
+## TODO: sort-based lookup is significantly slower given a large table
+## TODO: and small lookup vector.  This shouldn't be a problem since
+## TODO: interpolation (the reason for the table lookup in the first
+## TODO: place) usually involves subsampling of an existing table.  The
+## TODO: other use of interpolation, looking up values one at a time, is
+## TODO: unfortunately significantly slower for large tables.  
+## TODO:    sort is order O((lt+lx)*log(lt+lx)) 
+## TODO:    search is order O(lx*log(lt))
+## TODO: Clearly, search is asymptotically better than sort, but sort
+## TODO: is compiled and search is not.  Could support both, or recode
+## TODO: search in C++, or assume things are good enough as they stand.
+
+function idx = lookup (table, xi)
+  if (nargin == 2)
+    if (isempty (table))
+      idx = zeros (size (xi));
+    elseif (isvector (table))
+      [nr, nc] = size (xi);
+      lt = length (table);
+      if (table(1) > table(lt))
+	## decreasing table
+	[v, p] = sort ([xi(:); table(:)]);
+	idx(p) = cumsum (p > nr*nc);
+	idx = lt - idx(1:nr*nc);
+      else
+	## increasing table
+	[v, p] = sort ([table(:); xi(:) ]);
+	idx(p) = cumsum (p <= lt);
+	idx = idx(lt+1:lt+nr*nc);
+      endif
+      idx = reshape (idx, nr, nc);
+    else
+      error ("lookup: table must be a vector");
+    endif
+  else
+    print_usage ();
+  endif
+endfunction
+  
+%!assert (lookup(1:3, 0.5), 0)     # value before table
+%!assert (lookup(1:3, 3.5), 3)     # value after table error
+%!assert (lookup(1:3, 1.5), 1)     # value within table error
+%!assert (lookup(1:3, [3,2,1]), [3,2,1])
+%!assert (lookup([1:4]', [1.2, 3.5]'), [1, 3]');
+%!assert (lookup([1:4], [1.2, 3.5]'), [1, 3]');
+%!assert (lookup([1:4]', [1.2, 3.5]), [1, 3]);
+%!assert (lookup([1:4], [1.2, 3.5]), [1, 3]);
+%!assert (lookup(1:3, [3, 2, 1]), [3, 2, 1]);
+%!assert (lookup([3:-1:1], [3.5, 3, 1.2, 2.5, 2.5]), [0, 1, 2, 1, 1])
+%!assert (isempty(lookup([1:3], [])))
+%!assert (isempty(lookup([1:3]', [])))
+%!assert (lookup(1:3, [1, 2; 3, 0.5]), [1, 2; 3, 0]);

new file mode 100644
--- /dev/null
+++ b/scripts/general/nthroot.m
@@ -0,0 +1,71 @@
+## Copyright (C) 2004 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+##
+## Original version by Paul Kienzle distributed as free software in the
+## public domain.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} nthroot (@var{x}, @var{n})
+## 
+## Compute the nth root of @var{x}, returning real results for real 
+## components of @var{x}. For example
+##
+## @example
+## @group
+## nthroot (-1, 3)
+## @result{} -1
+## (-1) ^ (1 / 3)
+## @result{} 0.50000 - 0.86603i
+## @end group
+## @end example
+##
+## @end deftypefn
+
+function y = nthroot (x, m)
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+  
+  y = x.^(1./m);
+
+  if (isscalar (x))
+    x *= ones (size (m)); 
+  endif
+
+  if (isscalar (m))
+    m *= ones (size (x)); 
+  endif
+
+  idx = (mod (m, 2) == 1 & imag (x) == 0 & x < 0);
+
+  if (any (idx(:)))
+    y(idx) = -(-x(idx)).^(1./m(idx)); 
+  endif
+
+  ## If result is all real, make sure it looks real
+  if (all (imag (y) == 0))
+    y = real (y); 
+  endif
+
+endfunction
+
+%!assert(nthroot(-1,[3,-3]), [-1,-1],eps);
+%!assert(nthroot([-1,1],[3.1,-3]), [-1,1].^(1./[3.1,-3]));
+%!assert(nthroot([-1+1i,-1-1i],3), [-1+1i,-1-1i].^(1/3));

new file mode 100644
--- /dev/null
+++ b/scripts/linear-algebra/rref.m
@@ -0,0 +1,84 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{r}, @var{k}] =} rref (@var{a}, @var{tol})
+##
+## Returns the reduced row echelon form of @var{a}. @var{tol} defaults
+## to @code{eps * max (size (@var{a})) * norm (@var{a}, inf)}.
+##
+## Called with two return arguments, @var{k} returns the vector of
+## "bound variables", which are those columns on which elimination 
+## has been performed.
+##
+## @end deftypefn
+
+## Author: Paul Kienzle <pkienzle@users.sf.net>
+##         (based on a anonymous source from the public domain)
+
+function [A, k] = rref (A, tolerance)
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+  endif
+
+  if (ndims (A) > 2)
+    error ("rref: expecting matrix argument");
+  endif
+
+  [rows, cols] = size (A);
+
+  if (nargin < 2)
+    tolerance = eps * max (rows, cols) * norm (A, inf);
+  endif
+
+  used = zeros (1, cols);
+  r = 1;
+  for c = 1:cols
+    ## Find the pivot row
+    [m, pivot] = max (abs (A(r:rows,c)));
+    pivot = r + pivot - 1;
+
+    if (m <= tolerance)
+      ## Skip column c, making sure the approximately zero terms are
+      ## actually zero.
+      A (r:rows, c) = zeros (rows-r+1, 1);
+    else
+      ## keep track of bound variables
+      used (1, c) = 1;
+
+      ## Swap current row and pivot row
+      A ([pivot, r], c:cols) = A ([r, pivot], c:cols);
+
+      ## Normalize pivot row
+      A (r, c:cols) = A (r, c:cols) / A (r, c);
+
+      ## Eliminate the current column
+      ridx = [1:r-1, r+1:rows];
+      A (ridx, c:cols) = A (ridx, c:cols) - A (ridx, c) * A(r, c:cols);
+
+      ## Check if done
+      if (r++ == rows)
+	break;
+      endif
+    endif
+  endfor
+  k = find (used);
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/factor.m
@@ -0,0 +1,94 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{p} =} factor (@var{q})
+## @deftypefnx {Function File} {[@var{p}, @var{n}] =} factor (@var{q})
+##
+## Return prime factorization of @var{q}. That is @code{prod (@var{p})
+## == @var{q}}. If @code{@var{q} == 1}, returns 1. 
+##
+## With two output arguments, returns the uniques primes @var{p} and
+## their mulyiplicities. That is @code{prod (@var{p} .^ @var{n}) ==
+## @var{q}).
+## 
+## @end deftypefn
+
+## Author: Paul Kienzle
+
+## 2002-01-28 Paul Kienzle
+## * remove recursion; only check existing primes for multiplicity > 1
+## * return multiplicity as suggested by Dirk Laurie
+## * add error handling
+
+function [x, m] = factor (n)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  if (! isscalar (n) || n != fix (n))
+    error ("factor: n must be a scalar integer");
+  endif
+
+  ## special case of no primes less than sqrt(n)
+  if (n < 4)
+    x = n;
+    m = 1;
+    return;
+  endif 
+
+  x = [];
+  ## There is at most one prime greater than sqrt(n), and if it exists,
+  ## it has multiplicity 1, so no need to consider any factors greater
+  ## than sqrt(n) directly. [If there were two factors p1, p2 > sqrt(n),
+  ## then n >= p1*p2 > sqrt(n)*sqrt(n) == n. Contradiction.]
+  p = primes (sqrt (n));
+  while (n > 1)
+    ## find prime factors in remaining n
+    q = n ./ p;
+    p = p (q == fix (q));
+    if (isempty (p))
+      p = n;  # can't be reduced further, so n must itself be a prime.
+    endif
+    x = [x, p];
+    ## reduce n
+    n = n / prod (p);
+  endwhile
+  x = sort (x);
+
+  ## determine muliplicity
+  if (nargout > 1)
+    idx = find ([0, x] != [x, 0]);
+    x = x(idx(1:length(idx)-1));
+    m = diff (idx);
+  endif
+
+endfunction
+
+## test:
+##   assert(factor(1),1);
+##   for i=2:20
+##      p = factor(i);
+##      assert(prod(p),i);
+##      assert(all(isprime(p)));
+##      [p,n] = factor(i);
+##      assert(prod(p.^n),i);
+##      assert(all([0,p]!=[p,0]));
+##   end

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/isprime.m
@@ -0,0 +1,53 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} isprime (@var{n})
+##
+## Return true if @var{n} is a prime number, false otherwise.
+##
+## Something like the following is much faster if you need to test a lot
+## of small numbers:
+##
+## @example
+##    @var{t} = ismember (@var{n}, primes (max (@var{n} (:))));
+## @end example
+##
+## If max(n) is very large, then you should be using special purpose 
+## factorization code.
+##
+## @seealso{primes, factor, gcd, lcm}
+## @end deftypefn
+
+function t = isprime (n)
+  if (! isscalar (n))
+    nel = numel (n);
+    t = n;
+    for i = 1:nel
+      t(i) = isprime (t(i));
+    endfor
+  elseif (n != fix (n) || n < 2)
+    t = 0;
+  elseif (n < 9)
+    t = all (n != [4, 6, 8]);
+  else
+    q = n./[2, 3:2:sqrt(n)];
+    t = all (q != fix (q));
+  endif
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/legendre.m
@@ -0,0 +1,133 @@
+## Copyright (C) 2000  Kai Habel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{L} =} legendre (@var{n}, @var{X})
+##
+## Legendre Function of degree n and order m
+## where all values for m = 0..@var{n} are returned.
+## @var{n} must be a scalar in the range [0..255].
+## The return value has one dimension more than @var{x}.
+##
+## @example
+## The Legendre Function of degree n and order m
+##
+## @group
+##  m        m       2  m/2   d^m
+## P(x) = (-1) * (1-x  )    * ----  P (x)
+##  n                         dx^m   n
+## @end group
+##
+## with:
+## Legendre Polynom of degree n
+##
+## @group
+##           1     d^n   2    n
+## P (x) = ------ [----(x - 1)  ] 
+##  n      2^n n!  dx^n
+## @end group
+##
+## legendre(3,[-1.0 -0.9 -0.8]) returns the matrix
+##
+## @group
+##  x  |   -1.0   |   -0.9   |  -0.8
+## ------------------------------------
+## m=0 | -1.00000 | -0.47250 | -0.08000
+## m=1 |  0.00000 | -1.99420 | -1.98000
+## m=2 |  0.00000 | -2.56500 | -4.32000
+## m=3 |  0.00000 | -1.24229 | -3.24000 
+## @end group
+## @end example
+## @end deftypefn
+
+## FIXME Add Schmidt semi-normalized and fully normalized legendre functions
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+
+function L = legendre (n, x)
+
+  warning ("legendre is unstable for higher orders");
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+
+    if (! isscalar (n) || n < 0 || n > 255 || n != fix (n))
+      error ("n must be a integer between 0 and 255]");
+    endif
+
+    if (! isvector (x) || any (x < -1 || x > 1))
+      error ("x must be vector in range -1 <= x <= 1");
+    endif
+
+    if (n == 0)
+      L = ones (size (x));
+    elseif (n == 1)
+      L = [x; -sqrt(1 - x .^ 2)];
+    else
+      i = 0:n;
+      a = (-1) .^ i .* bincoeff (n, i);
+      p = [a; zeros(size (a))];
+      p = p(:);
+      p(length (p)) = [];
+      #p contains the polynom (x^2-1)^n
+
+      #now create a vector with 1/(2.^n*n!)*(d/dx).^n
+      d = [((n + rem(n, 2)):-1:(rem (n, 2) + 1)); 2 * ones(fix (n / 2), n)];
+      d = cumsum (d);
+      d = fliplr (prod (d'));
+      d = [d; zeros(1, length (d))];
+      d = d(1:n + 1) ./ (2 ^ n *prod (1:n));
+
+      Lp = d' .* p(1:length (d));
+      #Lp contains the Legendre Polynom of degree n
+
+      # now create a polynom matrix with d/dx^m for m=0..n
+      d2 = flipud (triu (ones (n)));
+      d2 = cumsum (d2);
+      d2 = fliplr (cumprod (flipud (d2)));
+      d3 = fliplr (triu (ones (n + 1)));
+      d3(2:n + 1, 1:n) = d2;
+
+      # multiply for each m (d/dx)^m with Lp(n,x)
+      # and evaluate at x
+      Y = zeros(n + 1, length (x));
+      [dr, dc] = size (d3);
+      for m = 0:dr - 1
+        Y(m + 1, :) = polyval (d3(m + 1, 1:(dc - m)) .* Lp(1:(dc - m))', x)(:)';
+      endfor
+
+      # calculate (-1)^m*(1-x^2)^(m/2)	for m=0..n at x
+      # and multiply with (d/dx)^m(Pnx)
+      l = length (x);
+      X = kron ((1 - x(:) .^ 2)', ones (n + 1, 1));
+      M = kron ((0:n)', ones (1, l));
+      L = X .^ (M / 2) .* (-1) .^ M .* Y;
+    endif
+endfunction
+
+%!test
+%! result=legendre(3,[-1.0 -0.9 -0.8]);
+%! expected = [
+%!    -1.00000  -0.47250  -0.08000
+%!     0.00000  -1.99420  -1.98000
+%!     0.00000  -2.56500  -4.32000
+%!     0.00000  -1.24229  -3.24000
+%! ];
+%! assert(result,expected,1e-5);

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/nchoosek.m
@@ -0,0 +1,78 @@
+## Copyright (C) 2001 Rolf Fabian and Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{c} =} nchoosek (@var{n}, @var{k})
+##
+## Compute the binomial coeeficient or all combinations of @var{n}.
+## If @var{n} is a scalar then, calculate the binomial coefficient
+## of @var{n} and @var{k}, defined as
+##
+## @iftex
+## @tex
+## $$
+##  {n \choose k} = {n (n-1) (n-2) \cdots (n-k+1) \over k!}
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+##
+## @example
+## @group
+##  /   \
+##  | n |    n (n-1) (n-2) ... (n-k+1)
+##  |   |  = -------------------------
+##  | k |               k!
+##  \   /
+## @end group
+## @end example
+## @end ifinfo
+##
+## If @var{n} is a vector generate all combinations of the elements
+## of @var{n}, taken @var{k} at a time, one row per combination. The 
+## resulting @var{c} has size @code{[nchoosek (length (@var{n}), 
+## @var{k}), @var{k}]}.
+##
+## @end deftypefn
+
+##AUTHORS Rolf Fabian  <fabian@tu-cottbus.de>
+##        Paul Kienzle <pkienzle@users.sf.net>
+
+## XXX FIXME XXX This function is identical to bincoeff for scalar
+## values, and so should probably be combined with bincoeff.
+
+function A = nchoosek (v, k)
+
+  n = length (v);
+
+  if (n == 1)
+     A = round (exp (sum (log (k+1:v)) - sum (log (2:v-k))));
+  elseif (k == 0)
+    A = [];
+  elseif (k == 1)
+    A = v(:);
+  elseif (k == n)
+     A = v(:).';
+  else
+    m = round (exp (sum (log (k:n-1)) - sum (log (2:n-k))));
+    A = [v(1)*ones(m,1), nchoosek(v(2:n),k-1);
+	 nchoosek(v(2:n),k)];
+  endif
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/perms.m
@@ -0,0 +1,47 @@
+## Copyright (C) 2001 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} perms (@var{v})
+##
+## Generate all permutations of @var{v}, one row per permutation. The
+## result has size @code{factorial (@var{n}) * @var{n}}, where @var{n}
+## is the length of @var{v}.
+##
+## @end deftypefn
+
+function A = perms (v)
+  v = v(:);
+  n = length (v);
+  if (n == 1)
+    A = v;
+  else
+    B = perms (v(1:n-1));
+    Bidx = 1:size (B, 1);
+    A = v(n) * ones (prod (2:n), n);
+    A(Bidx,1:n-1) = B;
+    k = size (B, 1);
+    for i = n-1:-1:2
+      A(k+Bidx,1:i-1) = B(Bidx,1:i-1);
+      A(k+Bidx,i+1:n) = B(Bidx,i:n-1);
+      k = k + size (B, 1);
+    endfor
+    A(k+Bidx,2:n) = B;
+  endif
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/specfun/primes.m
@@ -0,0 +1,88 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} primes (@var{n})
+##
+## Return all primes up to @var{n}.  
+##
+## Note that if you need a specific number of primes, you can use the
+## fact the distance from one prime to the next is on average
+## proportional to the logarithm of the prime.  Integrating, you find
+## that there are about @math{k} primes less than @math{k \log ( 5 k )}.
+##
+## The algorithm used is called the Sieve of Erastothenes.
+## @end deftypefn
+
+## Author: Paul Kienzle
+## Author: Francesco Potort́
+## Author: Dirk Laurie
+
+function x = primes (p)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  if (! isscalar (p))
+    error ("primes: n must be a scalar");
+  endif
+
+  if (p > 100000)
+    ## optimization: 1/6 less memory, and much faster (asymptotically)
+    ## 100000 happens to be the cross-over point for Paul's machine;
+    ## below this the more direct code below is faster.  At the limit
+    ## of memory in Paul's machine, this saves .7 seconds out of 7 for
+    ## p=3e6.  Hardly worthwhile, but Dirk reports better numbers.
+    lenm = floor ((p+1)/6);       # length of the 6n-1 sieve
+    lenp = floor ((p-1)/6);       # length of the 6n+1 sieve
+    sievem = ones (1, lenm);      # assume every number of form 6n-1 is prime
+    sievep = ones (1, lenp);      # assume every number of form 6n+1 is prime
+
+    for i = 1:(sqrt(p)+1)/6         # check up to sqrt(p)
+      if (sievem(i))              # if i is prime, eliminate multiples of i
+        sievem(7*i-1:6*i-1:lenm) = 0;
+        sievep(5*i-1:6*i-1:lenp) = 0;
+      endif                       # if i is prime, eliminate multiples of i
+      if (sievep(i))
+        sievep(7*i+1:6*i+1:lenp) = 0;
+        sievem(5*i+1:6*i+1:lenm) = 0;
+      endif
+    endfor
+    x = sort([2, 3, 6*find(sievem)-1, 6*find(sievep)+1]);
+  elseif (p > 352) # nothing magical about 352; just has to be greater than 2
+    len = floor ((p-1)/2);        # length of the sieve
+    sieve = ones (1, len);        # assume every odd number is prime
+    for i = 1:(sqrt(p)-1)/2       # check up to sqrt(p)
+      if (sieve(i))               # if i is prime, eliminate multiples of i
+        sieve(3*i+1:2*i+1:len) = 0; # do it
+      endif
+    endfor
+    x = [2, 1+2*find(sieve)];     # primes remaining after sieve
+  else
+    a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, ...
+	 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, ...
+	 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, ...
+	 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, ...
+	 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, ...
+	 293, 307, 311, 313, 317, 331, 337, 347, 349];
+    x = a(a <= p);
+  endif
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/special-matrix/hadamard.m
@@ -0,0 +1,173 @@
+## Copyright (C) 2004 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+##
+## Original version by Paul Kienzle distributed as free software in the
+## public domain.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} hadamard (@var{n})
+##
+## Construct a Hadamard matrix @var{Hn} of size @var{n}-by-@var{n}. The 
+## size @var{n} must be of the form @code{2 ^ @var{k} * @var{p}} where
+## @var{p} is one of 1, 12, 20 or 28. The returned matrix is normalized,
+## meaning @var{Hn (:, 1) == 1} and @var{H (1, :) == 1}.
+##
+## Some of the properties of Hadamard matrices are:
+##
+## @itemize @bullet
+## @item
+## @code{kron (@var{Hm}, @var{Hn})} is a Hadamard matrix of size 
+## @var{m}-by-@var{n}.
+## @item
+## @code{Hn * Hn' == @var{n) * eye (@var{n})}.
+## @item
+## The rows of @var{Hn} are orthogonal.
+## @item
+## @code{det (@var{A}) <= det (@var{Hn})} for all @var{A} with 
+## @code{abs (@var{A} (@var{i}, @var{j})) <= 1}.
+## @item
+## Multiply any row or column by -1 and still have a Hadamard matrix.
+## @end itemize
+##
+## @end deftypefn
+
+   
+## Reference [1] contains a list of Hadamard matrices up to n=256.
+## See code for h28 in hadamard.m for an example of how to extend
+## this function for additional p.
+##
+## References:
+## [1] A Library of Hadamard Matrices, N. J. A. Sloane 
+##     http://www.research.att.com/~njas/hadamard/
+
+function h = hadamard (n)
+  
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  ## find k if n = 2^k*p
+  k = 0;
+  while (n > 1 && floor (n/2) == n/2)
+    k++; 
+    n = n/2; 
+  endwhile
+  
+  ## find base hadamard
+  ## except for n=2^k, need a multiple of 4
+  if (n != 1)
+    k -= 2; 
+  endif
+
+  ## trigger error if not a multiple of 4
+  if (k < 0)
+    n =- 1;
+  endif
+
+  switch (n)
+    case 1
+      h = 1;
+    case 3
+      h = h12 ();
+    case 5
+      h = h20 ();
+    case 7
+      h = hnormalize (h28 ());
+    otherwise
+      error ("n must be 2^k*p, for p = 1, 12, 20 or 28");
+  endswitch
+
+  ## build H(2^k*n) from kron(H(2^k),H(n))
+  h2 = [1,1;1,-1];
+  while (true)
+    if (floor (k/2) != k/2)
+      h = kron (h2, h); 
+    endif
+    k = floor (k/2);
+    if (k == 0) 
+      break; 
+    endif
+    h2 = kron (h2, h2);
+  endwhile
+endfunction
+
+function h = h12 ()
+  tu=[-1,+1,-1,+1,+1,+1,-1,-1,-1,+1,-1];
+  tl=[-1,-1,+1,-1,-1,-1,+1,+1,+1,-1,+1];
+  ## note: assert(tu(2:end),tl(end:-1:2))
+  h = ones (12);
+  h(2:end,2:end) = toeplitz (tu, tl);
+endfunction
+
+
+function h = h20 ()
+  tu=[+1,-1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1,-1,-1,+1,+1,-1,-1];
+  tl=[+1,-1,-1,+1,+1,-1,-1,-1,-1,+1,-1,+1,-1,+1,+1,+1,+1,-1,-1];
+  ## note: assert(tu(2:end),tl(end:-1:2))
+  h = ones (20);
+  h(2:end,2:end) = fliplr (toeplitz (tu, tl));
+endfunction
+
+function h = hnormalize (h)
+  ## Make sure each row and column starts with +1
+  h(h(:,1)==-1,:) *= -1;
+  h(:,h(1,:)==-1) *= -1;
+endfunction
+
+function h = h28 ()
+  ## Williamson matrix construction from
+  ## http://www.research.att.com/~njas/hadamard/had.28.will.txt
+  s = ['+------++----++-+--+-+--++--';
+       '-+-----+++-----+-+--+-+--++-';
+       '--+-----+++---+-+-+----+--++';
+       '---+-----+++---+-+-+-+--+--+';
+       '----+-----+++---+-+-+++--+--';
+       '-----+-----++++--+-+--++--+-';
+       '------++----++-+--+-+--++--+';
+       '--++++-+-------++--+++-+--+-';
+       '---++++-+-----+-++--+-+-+--+';
+       '+---+++--+----++-++--+-+-+--';
+       '++---++---+----++-++--+-+-+-';
+       '+++---+----+----++-++--+-+-+';
+       '++++--------+-+--++-++--+-+-';
+       '-++++--------+++--++--+--+-+';
+       '-+-++-++--++--+--------++++-';
+       '+-+-++--+--++--+--------++++';
+       '-+-+-++--+--++--+----+---+++';
+       '+-+-+-++--+--+---+---++---++';
+       '++-+-+-++--+------+--+++---+';
+       '-++-+-+-++--+------+-++++---';
+       '+-++-+---++--+------+-++++--';
+       '-++--++-+-++-+++----++------';
+       '+-++--++-+-++-+++-----+-----';
+       '++-++---+-+-++-+++-----+----';
+       '-++-++-+-+-+-+--+++-----+---';
+       '--++-++++-+-+----+++-----+--';
+       '+--++-+-++-+-+----+++-----+-';
+       '++--++-+-++-+-+----++------+'];
+  h = (s=='+');
+  h(!h) = -1;
+endfunction
+
+%!assert(hadamard(1),1)
+%!assert(hadamard(2),[1,1;1,-1])
+%!test
+%!  for n=[1,2,4,8,12,24,48,20,28,2^9]
+%!    h=hadamard(n); assert(norm(h*h'-n*eye(n)),0);
+%!  end

new file mode 100644
--- /dev/null
+++ b/scripts/special-matrix/magic.m
@@ -0,0 +1,87 @@
+## Copyright (C) 1999-2001 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} magic (@var{n})
+##
+## Create an @var{n}-by-@var{n} magic square. Note that @code{magic
+## (@var{2})} is undefined since there is no 2-by-2 magic square.
+##
+## @end deftypefn
+
+function A = magic(n)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  if (n != floor (n) || n < 0 || n == 2)
+    error ("magic: n must be an positive integer not equal to 2");
+  endif
+
+  if (n == 0)
+
+    A = [];
+
+  elseif (mod (n, 2) == 1)
+
+    shift = floor ((0:n*n-1)/n);
+    c = mod ([1:n*n] - shift + (n-3)/2, n);
+    r = mod ([n*n:-1:1] + 2*shift, n);
+    A (c*n+r+1) = 1:n*n;
+    A = reshape (A, n, n);
+
+  elseif (mod (n, 4) == 0)
+
+    A = reshape (1:n*n, n, n)';
+    I = [1:4:n, 4:4:n];
+    J = fliplr (I);
+    A(I,I) = A(J,J);
+    I = [2:4:n, 3:4:n];
+    J = fliplr (I);
+    A(I,I) = A(J,J);
+
+  elseif (mod (n, 4) == 2)
+
+    m = n/2;
+    A = magic (m);
+    A = [A, A+2*m*m; A+3*m*m, A+m*m];
+    k = (m-1)/2;
+    if (k>1)
+      I = 1:m;
+      J = [2:k, n-k+2:n];
+      A([I,I+m],J) = A([I+m,I],J);
+    endif
+    I = [1:k, k+2:m];
+    A([I,I+m],1) = A([I+m,I],1);
+    I = k + 1;
+    A([I,I+m],I) = A([I+m,I],I);
+  
+  endif
+
+endfunction
+
+%!test
+%! for i=3:30
+%!   A=magic(i);
+%!   assert(norm(diff([sum(diag(A)),sum(diag(flipud(A))),sum(A),sum(A')])),0)
+%! endfor
+%!assert(isempty(magic(0)));
+%!assert(magic(1),1);
+%!error magic(2)

new file mode 100644
--- /dev/null
+++ b/scripts/special-matrix/pascal.m
@@ -0,0 +1,73 @@
+## Copyright (C) 1999 Peter Ekberg
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} pascal (@var{n}, @var{t})
+##
+## Return the Pascal matrix of order @var{n} if @code{@var{t} = 0}.
+## @var{t} defaults to 0. Return lower triangular Cholesky factor of 
+## the Pascal matrix if @code{@var{t} = 1}. This matrix is its own
+## inverse, that is @code{pascal (@var{n}, 1) ^ 2 == eye (@var{n})}.
+## If @code{@var{t} = 2}, return a transposed and  permuted version of
+## @code{pascal (@var{n}, 1)}, which is the cube-root of the identity
+## matrix. That is @code{pascal (@var{n}, 2) ^ 3 == eye (@var{n})}.
+##
+## @seealso{hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz
+##           hadamard, wilkinson, compan, rosser}
+## @end deftypefn
+
+## Author: Peter Ekberg
+##         (peda)
+
+function retval = pascal (n, t)
+
+  if (nargin > 2) || (nargin == 0)
+    print_usage ();
+  endif
+
+  if (nargin == 1)
+    t = 0;
+  endif
+
+  if (! is_scalar (n) || ! is_scalar (t))
+    error ("pascal: expecting scalar arguments, found something else");
+  endif
+
+  retval = diag ((-1).^[0:n-1]);
+  retval(:,1) = ones (n, 1);
+
+  for j = 2:n-1
+    for i = j+1:n
+      retval(i,j) = retval(i-1,j) - retval(i-1,j-1);
+    endfor
+  endfor
+
+  if (t == 0)
+    retval = retval*retval';
+  elseif (t == 2)
+    retval = retval';
+    retval = retval(n:-1:1,:);
+    retval(:,n) = -retval(:,n);
+    retval(n,:) = -retval(n,:);
+    if (rem(n,2) != 1)
+      retval = -retval;
+    endif
+  endif
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/special-matrix/rosser.m
@@ -0,0 +1,47 @@
+## Copyright (C) 1999 Peter Ekberg
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} rosser ()
+##
+## Returns the Rosser matrix. This is a difficult test case used to test
+## eigenvalue algorithms.
+##
+## @seealso{hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz
+##           hadamard, wilkinson, compan, pascal}
+## @end deftypefn
+
+## Author: Peter Ekberg
+##         (peda)
+
+function retval = rosser ()
+
+  if (nargin != 0)
+    print_usage ();
+  endif
+
+  retval = [611,   196,  -192,   407,    -8,   -52,   -49,    29;
+            196,   899,   113,  -192,   -71,   -43,    -8,   -44;
+           -192,   113,   899,   196,    61,    49,     8,    52;
+            407,  -192,   196,   611,     8,    44,    59,   -23;
+             -8,   -71,    61,     8,   411,  -599,   208,   208;
+            -52,   -43,    49,    44,  -599,   411,   208,   208;
+            -49,    -8,     8,    59,   208,   208,    99,  -911;
+             29,   -44,    52,   -23,   208,   208,  -911,    99];
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/special-matrix/wilkinson.m
@@ -0,0 +1,47 @@
+## Copyright (C) 1999 Peter Ekberg
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} wilkinson (@var{n})
+##
+## Return the Wilkinson matrix of order @var{n}.
+##
+## @seealso{hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz
+##           hadamard, rosser, compan, pascal}
+## @end deftypefn
+
+## Author: Peter Ekberg
+##         (peda)
+
+function retval = wilkinson (n)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  nmax = length (n);
+  if (! (nmax == 1))
+    error ("wilkinson: expecting scalar argument, found something else");
+  endif
+
+  side = ones (n-1, 1);
+  center = abs (-(n-1)/2:(n-1)/2);
+  retval = diag (side, -1) + diag (center) + diag (side, 1);
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/scripts/strings/strtok.m
@@ -0,0 +1,141 @@
+## Copyright (C) 2000 Paul Kienzle
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{tok}, @var{rem}] =} strtok (@var{str}, @var{delim})
+## 
+## Find all characters up to but not including the first character which
+## is in the string delim.  If @var{rem} is requested, it contains the
+## remainder of the string, starting at the first deliminator. Leading
+## delimiters are ignored.  If @var{delim} is not specified, space is assumed.
+##
+## @end deftypefn
+
+## TODO: check what to do for a null delimiter
+
+function [tok, rem] = strtok (str, delim)
+
+  if (nargin<1 || nargin > 2)
+    print_usage ();
+  endif
+
+  if (nargin < 2 || isempty (delim))
+    delim = " "; 
+  endif
+
+  if (isempty (str))
+    tok = rem = "";
+  elseif (length (delim) > 3)
+    start = 1;
+    len = length (str);
+    while (start <= len)
+      if (all (str(start) != delim))
+	break; 
+      endif
+      start++;
+    endwhile
+    stop = start;
+    while (stop <= len)
+      if (any (str(stop) == delim))
+	break; 
+      endif
+      stop++;
+    endwhile
+    tok = str(start:stop-1);
+    rem = str(stop:len);
+  else
+    if (length (delim) == 1)
+      idx = find (str == delim);
+    elseif (length (delim) == 2)
+      idx = find (str == delim(1) | str == delim(2));
+    else
+      idx = find (str == delim(1) | str == delim(2) | str == delim(3));
+    endif
+    if (isempty (idx))
+      tok = str;
+      rem = "";
+    else
+      ## Find first non-leading delimiter.
+      skip = find (idx(:)' != 1:length(idx));
+      if (isempty (skip))
+      	tok = str(idx(length(idx))+1:length(str));
+      	rem = "";
+      else
+      	tok = str(skip(1):idx(skip(1))-1);
+      	rem = str(idx(skip(1)):length(str));
+      endif
+    endif
+  endif
+
+endfunction
+
+%!demo
+%! strtok("this is the life")
+%! % split at the first space, returning "this"
+
+%!demo
+%! s = "14*27+31"
+%! while 1
+%!   [t,s] = strtok(s, "+-*/");
+%!   printf("<%s>", t);
+%!   if isempty(s), break; endif
+%!   printf("<%s>", s(1));
+%! endwhile
+%! printf("\n");
+%! % ----------------------------------------------------
+%! % Demonstrates processing of an entire string split on
+%! % a variety of delimiters. Tokens and delimiters are 
+%! % printed one after another in angle brackets.  The
+%! % string is:
+
+%!# test the tokens for all cases
+%!assert(strtok(""), "");             # no string
+%!assert(strtok("this"), "this");     # no delimiter in string
+%!assert(strtok("this "), "this");    # delimiter at end
+%!assert(strtok("this is"), "this");  # delimiter in middle
+%!assert(strtok(" this"), "this");    # delimiter at start
+%!assert(strtok(" this "), "this");   # delimiter at start and end
+%!assert(strtok(" "), ""(1:0));            # delimiter only
+
+%!# test the remainder for all cases
+%!test [t,r] = strtok(""); assert(r, "");
+%!test [t,r] = strtok("this"); assert(r, "");
+%!test [t,r] = strtok("this "); assert(r, " ");
+%!test [t,r] = strtok("this is"); assert(r, " is");
+%!test [t,r] = strtok(" this"); assert(r, "");
+%!test [t,r] = strtok(" this "); assert(r, " ");
+%!test [t,r] = strtok(" "); assert(r, "");
+
+%!# simple check with 2 and 3 delimeters
+%!assert(strtok("this is", "i "), "th");
+%!assert(strtok("this is", "ij "), "th");
+
+%!# test all cases for 4 delimiters since a different 
+%!# algorithm is used when more than 3 delimiters
+%!assert(strtok("","jkl "), "");
+%!assert(strtok("this","jkl "), "this");
+%!assert(strtok("this ","jkl "), "this");
+%!assert(strtok("this is","jkl "), "this");
+%!assert(strtok(" this","jkl "), "this");
+%!assert(strtok(" this ","jkl "), "this");
+%!assert(strtok(" ","jkl "), ""(1:0));
+
+%!# test 'bad' string orientations
+%!assert(strtok(" this "'), "this"');   # delimiter at start and end
+%!assert(strtok(" this "',"jkl "), "this"');