octave_packages/m/special-matrix/pascal.m

   1 ## Copyright (C) 1999-2012 Peter Ekberg
   2 ## Copyright (C) 2009 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {} pascal (@var{n})
  22 ## @deftypefnx {Function File} {} pascal (@var{n}, @var{t})
  23 ## Return the Pascal matrix of order @var{n} if @code{@var{t} = 0}.  @var{t}
  24 ## defaults to 0.  Return the pseudo-lower triangular Cholesky@tie{}factor of
  25 ## the Pascal matrix if @code{@var{t} = 1} (The sign of some columns may be
  26 ## negative).  This matrix is its own inverse, that is @code{pascal (@var{n},
  27 ## 1) ^ 2 == eye (@var{n})}.  If @code{@var{t} = -1}, return the true
  28 ## Cholesky@tie{}factor with strictly positive values on the diagonal.  If
  29 ## @code{@var{t} = 2}, return a transposed and permuted version of @code{pascal
  30 ## (@var{n}, 1)}, which is the cube root of the identity matrix.  That is,
  31 ## @code{pascal (@var{n}, 2) ^ 3 == eye (@var{n})}.
  32 ##
  33 ## @seealso{chol}
  34 ## @end deftypefn
  35
  36 ## Author: Peter Ekberg
  37 ##         (peda)
  38
  39 function retval = pascal (n, t = 0)
  40
  41   if (nargin < 1 || nargin > 2)
  42     print_usage ();
  43   elseif (! (isscalar (n) && isscalar (t)))
  44     error ("pascal: N and T must be scalars");
  45   elseif (! any (t == [-1, 0, 1, 2]))
  46     error ("pascal: expecting T to be -1, 0, 1, or 2, found %d", t);
  47   endif
  48
  49   retval = zeros (n);
  50   if (n > 0)
  51     retval(:,1) = 1;
  52   endif
  53
  54   if (t == -1)
  55     for j = 2:n
  56       retval(j:n,j) = cumsum (retval(j-1:n-1,j-1));
  57     endfor
  58   else
  59     for j = 2:n
  60       retval(j:n,j) = -cumsum (retval(j-1:n-1,j-1));
  61     endfor
  62   endif
  63
  64   if (t == 0)
  65     retval = retval*retval';
  66   elseif (t == 2)
  67     retval = rot90 (retval, 3);
  68     if (rem (n,2) != 1)
  69       retval *= -1;
  70     endif
  71   endif
  72
  73 endfunction
  74
  75
  76 %!assert (pascal (3,-1), [1,0,0;1,1,0;1,2,1])
  77 %!assert (pascal (3,0), [1,1,1;1,2,3;1,3,6])
  78 %!assert (pascal (3,0), pascal (3))
  79 %!assert (pascal (3,1), [1,0,0;1,-1,0;1,-2,1])
  80 %!assert (pascal (3,2), [1,1,1;-2,-1,0;1,0,0])
  81 %!assert (pascal (0,2), [])
  82
  83 %% Test input validation
  84 %!error pascal ()
  85 %!error pascal (1,2,3)
  86 %!error <N and T must be scalars> pascal ([1 2])
  87 %!error <N and T must be scalars> pascal (1, [1 2])
  88 %!error <expecting T to be> pascal (3,-2)
  89 %!error <expecting T to be> pascal (3,4)
  90