Mercurial > hg > octave-lyh

comparison scripts/special-matrix/pascal.m @ 9094:2e35cfcf6a6a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix, optimize & extend pascal
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 06 Apr 2009 11:16:20 +0200
parents 1bf0ce0930be
children 5235caf89e12
comparison
equal deleted inserted replaced
9093:8be05554bbd0 9094:2e35cfcf6a6a
1 ## Copyright (C) 1999, 2006, 2007 Peter Ekberg 1 ## Copyright (C) 1999, 2006, 2007 Peter Ekberg
2 ## Copyright (C) 2009 VZLU Prague
2 ## 3 ##
3 ## This file is part of Octave. 4 ## This file is part of Octave.
4 ## 5 ##
5 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by 7 ## under the terms of the GNU General Public License as published by
21 ## 22 ##
22 ## Return the Pascal matrix of order @var{n} if @code{@var{t} = 0}. 23 ## Return the Pascal matrix of order @var{n} if @code{@var{t} = 0}.
23 ## @var{t} defaults to 0. Return lower triangular Cholesky factor of 24 ## @var{t} defaults to 0. Return lower triangular Cholesky factor of
24 ## the Pascal matrix if @code{@var{t} = 1}. This matrix is its own 25 ## the Pascal matrix if @code{@var{t} = 1}. This matrix is its own
25 ## inverse, that is @code{pascal (@var{n}, 1) ^ 2 == eye (@var{n})}. 26 ## inverse, that is @code{pascal (@var{n}, 1) ^ 2 == eye (@var{n})}.
27 ## If @code{@var{t} = -1}, return its absolute value. This is the
28 ## standard pascal triangle as a lower-triangular matrix.
26 ## If @code{@var{t} = 2}, return a transposed and permuted version of 29 ## If @code{@var{t} = 2}, return a transposed and permuted version of
27 ## @code{pascal (@var{n}, 1)}, which is the cube-root of the identity 30 ## @code{pascal (@var{n}, 1)}, which is the cube-root of the identity
28 ## matrix. That is @code{pascal (@var{n}, 2) ^ 3 == eye (@var{n})}. 31 ## matrix. That is @code{pascal (@var{n}, 2) ^ 3 == eye (@var{n})}.
29 ## 32 ##
30 ## @seealso{hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz 33 ## @seealso{hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz
42 45
43 if (nargin == 1) 46 if (nargin == 1)
44 t = 0; 47 t = 0;
45 endif 48 endif
46 49
47 if (! is_scalar (n) || ! is_scalar (t)) 50 if (! isscalar (n) || ! isscalar (t))
48 error ("pascal: expecting scalar arguments, found something else"); 51 error ("pascal: expecting scalar arguments, found something else");
49 endif 52 endif
50 53
51 retval = diag ((-1).^[0:n-1]); 54 retval = zeros (n);
52 retval(:,1) = ones (n, 1); 55 retval(:,1) = 1;
53 56
54 for j = 2:n-1 57 if (t == -1)
55 for i = j+1:n 58 for j = 2:n-1
56 retval(i,j) = retval(i-1,j) - retval(i-1,j-1); 59 retval(j:n,j) = cumsum (retval (j-1:n-1,j-1));
57 endfor 60 endfor
58 endfor 61 else
62 for j = 2:n-1
63 retval(j:n,j) = -cumsum (retval (j-1:n-1,j-1));
64 endfor
65 endif
59 66
60 if (t == 0) 67 if (t == 0)
61 retval = retval*retval'; 68 retval = retval*retval';
62 elseif (t == 2) 69 elseif (t == 2)
63 retval = retval'; 70 retval = retval';