# HG changeset patch # User jwe # Date 1148279114 0 # Node ID 1fe78adb91bc2ca891e8b9b620a5672b31474da4 # Parent 6c6ff9b82577f3a7e2ec9785fdb2854f23e7cda1 [project @ 2006-05-22 06:25:14 by jwe] diff --git a/libcruft/ChangeLog b/libcruft/ChangeLog --- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,7 @@ +2006-05-22 John W. Eaton + + * lapack/dlantr.f, lapack/zlantr.f: New files. + 2006-05-03 David Bateman * lapack/dpocon.f, lapack/zpocon.f, lapack/dpotrs.f, diff --git a/scripts/ChangeLog b/scripts/ChangeLog --- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,19 @@ +2006-05-22 John W. Eaton + + * scripts/general/lookup.m: New file from Octave Forge. + +2006-05-22 David Bateman + + * 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 * polynomial/unmkpp.m, polynomial/mkpp.m, polynomial/spline.m, diff --git a/scripts/general/lookup.m b/scripts/general/lookup.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 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]); diff --git a/scripts/general/nthroot.m b/scripts/general/nthroot.m 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)); diff --git a/scripts/linear-algebra/rref.m b/scripts/linear-algebra/rref.m 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 +## (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 diff --git a/scripts/specfun/factor.m b/scripts/specfun/factor.m 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 diff --git a/scripts/specfun/isprime.m b/scripts/specfun/isprime.m 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 diff --git a/scripts/specfun/legendre.m b/scripts/specfun/legendre.m 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 + +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); diff --git a/scripts/specfun/nchoosek.m b/scripts/specfun/nchoosek.m 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 +## Paul Kienzle + +## 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 diff --git a/scripts/specfun/perms.m b/scripts/specfun/perms.m 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 diff --git a/scripts/specfun/primes.m b/scripts/specfun/primes.m 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 diff --git a/scripts/special-matrix/hadamard.m b/scripts/special-matrix/hadamard.m 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 diff --git a/scripts/special-matrix/magic.m b/scripts/special-matrix/magic.m 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) diff --git a/scripts/special-matrix/pascal.m b/scripts/special-matrix/pascal.m 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 diff --git a/scripts/special-matrix/rosser.m b/scripts/special-matrix/rosser.m 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 diff --git a/scripts/special-matrix/wilkinson.m b/scripts/special-matrix/wilkinson.m 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 diff --git a/scripts/strings/strtok.m b/scripts/strings/strtok.m 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"');