Re: problem with consistency of root finding and sorting

[David Bateman]

Ben Abbott wrote:
> On Oct 9, 2007, at 8:23 PM, David Bateman wrote:
> 
>> Ben Abbott wrote:
>>> I have Octave 2.9.14 on Mac OS X (intel & PPC) installed via Fink and
>>> am concerned that my installations are giving different results than
>>> others.
>>>
>>> During the process of modifying residue.m so that it would function
>>> in a manner consistent with Matlab's I've become suspicious that my
>>> Octave installation  does not work the same as others.
>>>
>>> The difference is small, but is troublesome.
>>>
>>> Given the following short example,
>>>
>>>     r = sort(roots([1 0 18 0 81]))
>>>     r_err = r - round(r)
>>>
>>> I obtain the results below
>>>
>>> r = sort(roots([1 0 18 0 81]))
>>> r =
>>>
>>>     0.0000 - 3.0000i
>>>     0.0000 + 3.0000i
>>>    -0.0000 - 3.0000i
>>>    -0.0000 + 3.0000i
>>>
>>> r_err = r - round(r)
>>> r_err =
>>>
>>>    2.5814e-08 + 1.6920e-08i
>>>    2.5814e-08 - 1.6920e-08i
>>>    -2.5814e-08 - 1.6920e-08i
>>>    -2.5814e-08 + 1.6920e-08i
>>
>> Under Octave 2.9.14 compiled with g++ 4.0.1 on a Mandriva 2006  
>> machine I get
>>
>> octave:3> r = sort(roots([1 0 18 0 81]))
>> r =
>>
>>    0.0000 - 3.0000i
>>    0.0000 + 3.0000i
>>   -0.0000 - 3.0000i
>>   -0.0000 + 3.0000i
>>
>> octave:4> r_err = r - round(r)
>> r_err =
>>
>>   1.5973e-08 + 6.5742e-09i
>>   1.5973e-08 - 6.5742e-09i
>>   -1.5973e-08 - 6.5742e-09i
>>   -1.5973e-08 + 6.5742e-09i
>>
>> D.
> 
> The order of your roots is the same as mine! I've attached two files,  
> residue.m and mpoles.m. Perhaps the small differences in r_err are  
> due to slight differences in BLAS? <shrugs>
> 
> David, if you have the time, please download /scripts/polynomial/ 
> residue.m and /scripts/polynomial/mpoles.m from the cvs, and try  
> another test for me.
> 
> b = [1 0 1];
> a = [1 0 18 0 81];
> [r, p, k, e] = residue(b, a)
> 
> r =
> 
>    -0.00000 - 0.09259i
>     0.22222 - 0.08333i
>     0.00000 + 0.09259i
>     0.22222 + 0.08333i
> 
> p =
> 
>     0.0000 + 3.0000i
>     0.0000 + 3.0000i
>    -0.0000 - 3.0000i
>    -0.0000 - 3.0000i
> 
> k = [](0x0)
> e =
> 
>     1
>     2
>     1
>     2


For the CVS of a couple of hours ago, I get

octave:1> b = [1 0 1];
octave:2> a = [1 0 18 0 81];
octave:3> [r, p, k, e] = residue(b, a)
r =

  -0.00000 + 0.09259i
   0.22222 - 0.00000i
  -0.00000 - 0.09259i
   0.22222 + 0.00000i

p =

  -0.0000 - 3.0000i
  -0.0000 - 3.0000i
   0.0000 + 3.0000i
   0.0000 + 3.0000i

k = [](0x0)
e =

   1
   2
   1
   2

D.