allroots patch



Raymond Toy wrote:
> Andrej Vodopivec wrote:
>   
>> The polynomial which can't be solved is
>> x^6-24576*x^5+402653184*x^4-4947802415966*x^3+40532397764222976*x^2+9157742690304*x+2069067169.
>>
>> I tried bfallroots with fpprec=32 and is still didn't find any roots.
>> With fpprec 64 it only finds 2 roots. With fpprec 128 it finds all
>> roots.
>>   
>>     
> Interesting.  I wonder what makes the polynomial so difficult to solve. 
> The roots are pretty well separated, except for a pair near the origin.
>   
I vaguely remember reading somewhere that sometimes it's better to use
the complex Jenkins-Traub algorithm, even if the polynomial has real
coefficients.  And indeed, this is the case for this polynomial.  If we
multiply the polynomial by %i, then allroots, and bfallroots (with
fpprec = 16) instantly finds all 6 roots.  The complex conjugates are
exactly conjugates, but pretty close.

Maybe it is worthwhile to modify the algorithm slightly so that the
complex Jenkins-Traub is used for real polynomials, but  for each root
found, we remove the complex conjugate right away too.  That might be
tricky if there's a real root but roundoff produces a very small
imaginary component.

> And for me with fpprec = 64, I get an error.  Maxima prints out the
> message that only two roots are found and then an error about 0 not
> being a symbol.  Need to look into that.
>   
Found the problem.  Fix will be checked in shortly.  An intermediate
result was supposed to be initialized to a bigfloat zero, but the
bigfloat marker was not included.

Ray