changeset 9094:2e35cfcf6a6a

fix, optimize & extend pascal
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 06 Apr 2009 11:16:20 +0200
parents 8be05554bbd0
children b5a1a2d0be7e
files scripts/ChangeLog scripts/special-matrix/pascal.m
diffstat 2 files changed, 18 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Mon Apr 06 08:33:49 2009 +0200
+++ b/scripts/ChangeLog	Mon Apr 06 11:16:20 2009 +0200
@@ -1,3 +1,7 @@
+2009-04-06  Jaroslav Hajek  <highegg@gmail.com>
+
+	* special-matrix/pascal.m: Fix, optimize & extend.
+
 2009-04-06  Jaroslav Hajek  <highegg@gmail.com>
 
 	* linear-algebra/vech.m: Optimize.
--- a/scripts/special-matrix/pascal.m	Mon Apr 06 08:33:49 2009 +0200
+++ b/scripts/special-matrix/pascal.m	Mon Apr 06 11:16:20 2009 +0200
@@ -1,4 +1,5 @@
 ## Copyright (C) 1999, 2006, 2007 Peter Ekberg
+## Copyright (C) 2009 VZLU Prague
 ##
 ## This file is part of Octave.
 ##
@@ -23,6 +24,8 @@
 ## @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} = -1}, return its absolute value. This is the
+## standard pascal triangle as a lower-triangular matrix.
 ## 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})}.
@@ -44,18 +47,22 @@
     t = 0;
   endif
 
-  if (! is_scalar (n) || ! is_scalar (t))
+  if (! isscalar (n) || ! isscalar (t))
     error ("pascal: expecting scalar arguments, found something else");
   endif
 
-  retval = diag ((-1).^[0:n-1]);
-  retval(:,1) = ones (n, 1);
+  retval = zeros (n);
+  retval(:,1) = 1;
 
-  for j = 2:n-1
-    for i = j+1:n
-      retval(i,j) = retval(i-1,j) - retval(i-1,j-1);
+  if (t == -1)
+    for j = 2:n-1
+      retval(j:n,j) = cumsum (retval (j-1:n-1,j-1));
     endfor
-  endfor
+  else
+    for j = 2:n-1
+      retval(j:n,j) = -cumsum (retval (j-1:n-1,j-1));
+    endfor
+  endif
 
   if (t == 0)
     retval = retval*retval';