Re: Octave and Freemat

[David Bateman]

David Bateman wrote:
> Sebastien Loisel wrote:
>> Dear David,
>>
>> Thank you for your quick reply.
>>
>>
>>     That is an error, as the code on octave-forge is not part of
>>     Octave. The
>>     copyright strings were recently updated in octave-forge and I believe
>>
>>
>> OK, I see that ORTH 
>> http://velveeta.che.wisc.edu/cgi-bin/cvsweb.cgi/~checkout~/octave/scripts/linear-algebra/orth.m?rev=HEAD&content-type=text/plain
>> <http://velveeta.che.wisc.edu/cgi-bin/cvsweb.cgi/%7Echeckout%7E/octave/scripts/linear-algebra/orth.m?rev=HEAD&content-type=text/plain>
>>
>> is part of Octave, while ODE45
>> http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/extra/ode/inst/ode45.m?revision=HEAD&content-type=text/plain
>> <http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/extra/ode/inst/ode45.m?revision=HEAD&content-type=text/plain>
>> is not. I'll just leave the copyrights alone. For those files that I
>> lift out of Octave (like ORTH), I will put a note that the file was
>> copied into Freemat, and that Freemat's not part of Octave, and
>> otherwise leave the copyright alone. Does that sound like the right idea?
>>  
>>
>>     roots is an octave core function and not from Octave forge, so
>>     send the
>>     proposed patches to roots to address@hidden
>>     <mailto:address@hidden>. As for octave-forge
>>
>>
>> The correct implementation of roots.m is
>>
>> function z = roots(p)
>>   if(any(isnan(p) | isinf(p)))
> Octave deliberately doesn't support the short circuiting with the "|"
> operator, use "||" instead. See
> 
> http://www.gnu.org/software/octave/FAQ.html#MATLAB-compatibility
> 
> for a discussion of why.
> 
>>  
>>   while(any(isinf(p./p(1))))
>>      p=p(2:end);
>>   end
> Its bad form to use a loop here, but this is the cause of the issue you
> saw. It would be better to write that as
> 
> f = find (p ./ max(p));
> p = p (f(1):end);
> 
> assuming p is a vector of cause, so perhaps you should move your "vec"
> statement before this. I'll supply a patch to Octave's roots function
> based on this
> 
> D.
> 
> 


Ok attached is  what I would suggest as a modification to the Octave
roots function to ix this issue based on your contributed code.

D.

# HG changeset patch
# User Sebastien Loisel
# Date 1204671561 -3600
# Node ID 899bdac2abcca2bb160d25d6b8287263c3985a81
# 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,8 @@ 2008-03-04  Bill Denney  <address@hidden
+2008-03-04  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
@@ -91,7 +91,11 @@ function r = roots (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);
+  if (isempty (v))
+    f = v;
+  else
+    f = find (v ./ max(v));
+  endif
   m = max (size (f));
 
   if (m > 0 && n > 1)
@@ -120,3 +124,4 @@ endfunction
 
 %!error roots ([1, 2; 3, 4]);
 
+%!assert(roots ([1e-200, 1e200, 1]), -1e-200)