Ben Abbott wrote:
On Mar 4, 2008, at 8:27 PM, Ben Abbott wrote:
On Mar 4, 2008, at 6:46 PM, Sebastien Loisel wrote:
Dear David,
Thank you for your email.
On Tue, Mar 4, 2008 at 5:44 PM, David Bateman
<address@hidden> wrote:
f = find (p ./ max(p));
p = p (f(1):end);
Are you missing an abs maybe? Also, I hope you did your checking
for
nans and infs before you got there.
--
Sébastien Loisel
To respect Matlab an error should result when NaNs or Infs are
present.
The abs shouldn't be necessary. In fact, if NaNs and Infs have
already be handled, why not
f = find (p);
p = p(f(1):end);
n = numel (p)-1;
A = diag (ones (n-1, 1), -1);
A(1,:) = -p(2:n+1) ./ p(1);
z = eig (A);
Ben
ok, nix the simple solution.
I checked Matlab, they apparently remove have some special handling
of
trailing zeros.
p = [poly([3 3 3 3]), 0 0 0 0];
abs(roots(p)-[0; 0; 0; 0; 3; 3; 3; 3])
ans =
0
0
0
0
0.00051228
0.00051228
0.00051223
0.00051223
What should be included is the check for Infs and NaNs. Beyond that
we
might add some tests for consistency with Matlab.
How about the attached?
How about the attached instead that is a combination of what Sebastien
pointed out and your fixes..
D.
--
David Bateman
address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified
as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
# HG changeset patch
# User Sebastien Loisel
# Date 1204708216 -3600
# Node ID 1dea613279ff9b360b5df72ed3779ca7fcb9f166
# Parent 1ea3f5aa672f0bb1919a4355d38fee00690028c2
Apply a scaling factor to leading zero removal in roots.m
diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,12 @@ 2008-03-04 Bill Denney <address@hidden
+2008-03-05 Bill Denney <address@hidden>
+
+ * polynomial/roots.m: Modified to catch Infs and/or NaNs.
+
+2008-03-05 Sebastien Loisel <address@hidden>
+
+ * polynomial/roots.m: Apply a scaling fator to the removal of the
+ leading zeros.
+
2008-03-04 Bill Denney <address@hidden>
* geometry/rectint.m: New function.
diff --git a/scripts/polynomial/roots.m b/scripts/polynomial/roots.m
--- a/scripts/polynomial/roots.m
+++ b/scripts/polynomial/roots.m
@@ -83,16 +83,22 @@ function r = roots (v)
if (nargin != 1 || min (size (v)) > 1)
print_usage ();
+ elseif (any (isnan(v) | isinf(v)))
+ error ("roots: inputs must not contain Inf or NaN.")
endif
- n = length (v);
- v = reshape (v, 1, n);
+ n = numel (v);
+ v = v(:);
## 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 (isempty (v))
+ f = v;
+ else
+ f = find (v ./ max (abs (v)));
+ endif
+ m = numel (f);
if (m > 0 && n > 1)
v = v(f(1):f(m));
@@ -114,9 +120,24 @@ function r = roots (v)
endfunction
+%!test
+%! p = [poly([3 3 3 3]), 0 0 0 0];
+%! r = sort (roots (p));
+%! assert (r, [0; 0; 0; 0; 3; 3; 3; 3], 0.001)
+
%!assert(all (all (abs (roots ([1, -6, 11, -6]) - [3; 2; 1]) < sqrt
(eps))));
%!assert(isempty (roots ([])));
%!error roots ([1, 2; 3, 4]);
+
+%!assert(isempty (roots (1)));
+ %!error roots ([1, 2; 3, 4]);
+
+%!error roots ([1 Inf 1]);
+
+%!error roots ([1 NaN 1]);
+
+%!assert(roots ([1e-200, -1e200, 1]), 1e-200)
+%!assert(roots ([1e-200, -1e200 * 1i, 1]), -1e-200 * 1i)