I've used an Aberth's method implementation by D. A. Bini:
http://www.netlib.org/numeralgo/na10
Here's my Lisp draft: http://beshenov.ru/maxima/pzeros.lisp
I didn't test it carefully (for example, it works fine on sbcl but
fails on clisp [1]), but I guess it could be extended to the arbitrary
precision and it will be also useful or interesting for someone...
A short example:
(polzeros #(1.0 2.0 3.0 4.0)) =>
(#(#C(-0.07208523 -0.6383267) #C(-0.07208521 0.63832676)
#C(-0.6058296 0.0))
#(1.4296476e-6 1.3818666e-6 1.8927716e-6))
The first array represents the approximations to the complex roots and
the second array represents the error bounds.
[1] As for clisp vs. sbcl, clisp 2.44.1 gives (on my machine)
least-positive-single-float => 1.1754944E-38
most-positive-single-float => 3.4028235E38
(/ 1.0 most-positive-single-float) => floating point underflow
(actually, 1.1754944E-38 is a least-positive-normalized-single-float!)
while sbcl 1.0.18 gives
least-positive-single-float => 1.4012985e-45
most-positive-single-float => 3.4028235e38
(/ 1.0 most-positive-single-float) => 2.938736e-39
On Friday 12 December 2008 22:01:32 you wrote:
> Some time ago, I wrote jtroot3.mac (share/numeric) to compute roots
> of polynomials using Jenkins-Traub algorithm. It mostly works, but
> some simple polynomials fail (see
> https://sourceforge.net/tracker2/?func=detail&aid=2230778&group_id=
>4933&atid=104933).
>
> What I should have done is taken the existing allroots
> implementation and converted it to use bfloats. Well, here it is:
> http://common-lisp.net/~rtoy/bfcpoly.lisp.
>
> There are probably bugs, but the few tests I tried seem to work,
> and the above bug is gone.
>
> While doing this, I also found the undocumented polyfactor
> variable. When polyfactor is false, bfallroots and allroots return
> a list of the roots. When polyfactor is true, bfallroots and
> allroots return the result as a factored polynomial of the form
> c*(x-r1)*(x-r2)...
>
> Perhaps this can be placed in cpoly.lisp or maybe as a contrib in
> share/numeric.
>
> Also, I should probably also convert rpoly-sl to support bfloats,
> and then allroots and bfallroots will be functionally identical.
--
Boomtime, Aftermath 55 YOLD 3174
Alexey Beshenov http://beshenov.ru/