1 ## Copyright (C) 1999-2012 Peter Ekberg
2 ## Copyright (C) 2009 VZLU Prague
4 ## This file is part of Octave.
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.
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.
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/>.
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})}.
36 ## Author: Peter Ekberg
39 function retval = pascal (n, t = 0)
41 if (nargin < 1 || nargin > 2)
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);
56 retval(j:n,j) = cumsum (retval(j-1:n-1,j-1));
60 retval(j:n,j) = -cumsum (retval(j-1:n-1,j-1));
65 retval = retval*retval';
67 retval = rot90 (retval, 3);
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), [])
83 %% Test input validation
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)