Octave and Freemat

[Ben Abbott]

On Mar 5, 2008, at 4:12 AM, David Bateman wrote:

> 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
>>>> <David.Bateman at motorola.com> 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                                 
> David.Bateman at motorola.com
> 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  <bill at denney.ws
> +2008-03-05  Bill Denney  <bill at denney.ws>
> +
> +	* polynomial/roots.m: Modified to catch Infs and/or NaNs.
> +
> +2008-03-05  Sebastien Loisel  <loisel at temple.edu>
> +
> +	* polynomial/roots.m: Apply a scaling fator to the removal of the
> +	leading zeros.
> +
> 2008-03-04  Bill Denney  <bill at denney.ws>
>
> 	* 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)

Regarding this part,

> +  if (isempty (v))
> +    f = v;
> +  else
> +    f = find (v ./ max (abs (v)));
> +  endif
> +  m = numel (f);

It appears an error will result if v = 0. Which should return [].  
Perhaps something like that below?

   f = find (v);
   if (any (f))
     v = v ./ max (abs (v));
     f = find (v);
   endif
   m = numel (f);

The first test below tests this concern, and the second tests for a  
"special" features of Matlab's implementation, which apparently is  
consistent with the above since the check for Inf is done prior to  
"v ./ max (abs (v))".

%!assert(isempty (roots (0)));

%!assert(roots([realmin, realmax, realmax]), -1)

Ben