--- a/scripts/polynomial/roots.m	Wed Nov 01 00:48:17 1995 +0000
+++ b/scripts/polynomial/roots.m	Wed Nov 01 01:35:40 1995 +0000
@@ -23,39 +23,38 @@
 # For a vector v with n components, return the roots of the
 # polynomial v(1) * z^(n-1) + ... + v(n-1) * z + v(n).
   
-# Written by KH (Kurt.Hornik@neuro.tuwien.ac.at) on Dec 24, 1993.
+# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993
+# Updated by KH on Nov 25, 1994
+  
+  if (min (size (v)) > 1 || nargin != 1)
+    usage ("roots (v), where v is a vector");
+  endif
+  
+  n = length (v);
+  v = reshape (v, 1, n);  
 
-  [nr, nc] = size (v);
-  if (nr <= 1 && nc <= 1)
-    r = [];
-    return;
-  elseif (! ((nr == 1 && nc > 1) || (nc == 1 && nr > 1)))
-    usage ("roots (v), where v is a nonzero vector");
-  endif
-
-  n = nr + nc - 1;
-  v = reshape (v, 1, n);
-
-# If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
-# leading k zeros and n - k - l roots of the polynomial are zero.  
+  # If v = [ 0 ... 0 v(k+1) ... v(k+l) 0 ... 0 ], we can remove the
+  # leading k zeros and n - k - l roots of the polynomial are zero.  
 
   f = find (v);
   m = max (size (f));
-  if (m > 0)
-    v = v (f(1):f(m));
+  
+  if (m > 0 && n > 1)
+    v = v(f(1):f(m));
     l = max (size (v));
     if (l > 1)
       A = diag (ones (1, l-2), -1);
-      A (1, :) = -v (2:l) ./ v (1);
+      A(1,:) = -v(2:l) ./ v(1);
       r = eig (A);    
-      if (f (m) < n)
-	r = [r; (zeros(n - f(m), 1))];
+      if (f(m) < n)
+        tmp = zeros (n - f(m), 1);
+	r = [r; tmp];
       endif
     else
       r = zeros (n - f(m), 1);
     endif
   else
-    usage ("roots (v), where v is a nonzero vector");
+    r = [];
   endif
   
 endfunction