--- tril.m.orig	2009-10-22 17:36:30.000000000 +0200
+++ tril.m	2009-10-22 17:49:14.000000000 +0200
@@ -72,22 +72,66 @@
      endif 
     [nr, nc] = size (x);
   endif
-
   if (nargin == 1)
     k = 0;
   elseif (nargin == 2)
-    if ((k > 0 && k > nc) || (k < 0 && k < -nr))
+    if ((k > 0 && k > nc-1) || (k < 0 && k < -nr))
       error ("tril: requested diagonal out of range");
     endif
   else
     print_usage ();
   endif
-
-  retval = resize (resize (x, 0), nr, nc);
-  for j = 1 : min (nc, nr+k)
-    nr_limit = max (1, j-k);
-    retval (nr_limit:nr, j) = x (nr_limit:nr, j);
-  endfor
+  if (nr == nc)
+    if (k < 0)
+      retval = zeros(size(x));      
+      for j = 1 : nr+k
+	idx = (j-1)*nr+j-k:j*nr;
+	retval(idx) = x(idx);
+      endfor	
+    else
+% k >= 0
+      retval = x;
+      for j = nc : -1 : k+2
+	retval((j-1)*nr+1:(j-1)*nr-k-1+j) = 0;
+      endfor	
+    endif
+  elseif (nr > nc)
+    if (k < floor((1+nc-nr)/2))
+      retval = zeros(size(x));      
+      for j = 1 : min(nc,nr+k)
+	idx = (j-1)*nr+j-k:j*nr;
+	retval(idx) = x(idx);
+      endfor	
+    else
+% k >= floor((1+nc-nr)/2)
+      retval = x;
+      for j = nc : -1 : max(1,k+2)
+	retval((j-1)*nr+1:(j-1)*nr-k-1+j) = 0;
+      endfor	
+    endif
+  else
+% nc > nr
+    if (k < floor((1+nc-nr)/2))
+      retval = zeros(size(x));      
+      for j = 1 : k
+	idx = (j-1)*nr+1:j*nr;
+	retval(idx) = x(idx);
+      endfor	
+      for j = max(1,k+1) : nr+k
+	idx = (j-1)*nr+j-k:j*nr;
+	retval(idx) = x(idx);
+      endfor	
+    else
+% k >= floor((1+nc-nr)/2)
+      retval = x;
+      for j = nc : -1 : k+nr+2
+	retval((j-1)*nr+1:j*nr) = 0;
+      end	
+      for j = min(nc,k+1+nr) : -1 : k+2
+	retval((j-1)*nr+1:(j-1)*nr-k-1+j) = 0;
+      endfor	
+    endif
+  endif
 
 endfunction
 
@@ -109,4 +153,3 @@
 %!error tril ();
 
 %!error tril (1, 2, 3);
-