oldk1331 wrote:
>
> Seems there's no way in FriCAS to compute root of polynomials in
> DoubleFloat?
>
> If so, we should add it, it should be much fater than Float.
> I wonder if it can be easily added in numsolve.spad.
Well, finding/creating good numeric solver is not easy.
We currently have real solver based on Descartes rule
of signs. It works well up to degree of order 100.
The method should work resonably well for degrees
of order 2000, but we would need to spend significat
work on improving efficiency.
There is alternative solver in ComplexRootPackage based
on Shoenage modification to Lobachevsky method. Theoretically
is should have quite good asymptotic complexity, but
its operation depends on correct choice of precision:
too low and it will fail to converge, too big and
it will be slow. To make things more tricky optimal
accuracy of various steps may differ. Currently this
solver assumes that user has selected correct
accuracy so it is hard to use. Also, it misses
various optimizatins.
Coming back to double precision:
- it is tricky to correctly evaluate high degree polynomial
in double precision
- we need error estimates, theoretically possible (and
not very hard), but adds to needed code
- numeric method may fail to converge, in particular it
seems that method which converges when using high
accuracy may fail to converge in lower accuracy
- convergence seem to depend quite a lot on initial
approximation
For some time Jenkins-Traub method was considerd to be
standard. But the code for Jenkins-Traub is quite
long and later research showed that it behaves
badly for polynomials of high degree (say above 20).
I experimented a bit with few different methods and
they seem to work much worse than advertised. To
put it differently: it seems that they work for
authors only because of various undisclosed tricks
(alternatively, authors did inadequate testing
and only report positive results). Particulary
discouraging is fact that IIRC in one case
method started converging only after I bumped
accuracy to quite high value -- that suggests
that in double precision the method is divergent.
There is related theoretical question: given
appropriatly normalized polynomial of degree n
how many bits of precion are needed to solve it?
The best result seem to be due to V. Pan and says
that C*n*log(n) bits is enough (from the proof
one gets C slightly larger than 7). This result
may be optimal, because several intermedate
steps get unstable with lower accuracy.
To summarize: it should be possible to create
good solver for low degree polynomials, but the
effort is nontrivial. In the literature there
were reports about solving polynomials of degrees
of order 1000 using double precision arithmetic,
but it is not clear if such solver can be general
or can only handle some well-behaved polynomials.
--
Waldek Hebisch
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
