Re: Octave and Freemat

[David Bateman]

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)