Hope nobody minds, but I'm appending a small maxima program to compute
the roots of a polynomial with bfloat coefficients. It uses
Jenkins-Traub to compute the roots, much like allroots. But this
doesn't have all the refinements of allroots. The code is simple, and
could probably be done much better---I'm not a very good maxima
programmer.
The polynomial is given as a list of the coefficients, arranged in
descending powers. The code should work with rational coefficients
or float coefficients too, and would probably converge, provided you
set fpprec to something reasonable, because fpprec is used to
determine if the algorithm has converged.
A quick example of computing the roots of a polynomial whose roots are
1 through 20, using 32 digits.
(%i230) fpprec:32$
(%i231) p:expand(product((x-k),k,1,20))$
(%i232) bfp:map(lambda([n], bfloat(coeff(p,x,20-n))), makelist(k,k,0,20))$
(%i233) polyroots(bfp);
(%o233) [6.5481618109166019003903250047888b-33 %i
+ 9.9999999999999999999999999999999b-1,
1.9999999999999999999999999999791b0 - 9.1674265352832426605464550067043b-32
%i, 1.8797076257219421925826344719629b-31 %i
+ 3.0000000000000000000000000024117b0,
5.03361050265047962554710395074b-30 %i + 3.9999999999999999999999999381496b0,
5.0000000000000000000000007275936b0 - 1.2563509556040235166793456910545b-25
%i, 1.8840508535100320240561969880573b-24 %i
+ 5.9999999999999999999999955144086b0,
7.0000000000000000000000114365974b0 - 1.3187685895779509228236651585887b-23
%i, 8.9999999999999999999994656020482b0
- 1.7144064831208246850001620131232b-22 %i,
5.7143431224331995655745113480833b-23 %i
+ 8.0000000000000000000000415222514b0,
1.0999999999999999999991786043188b1 - 6.2990391246214643613810666070531b-22
%i, 3.7725039187392241115037372420251b-22 %i
+ 1.0000000000000000000002628892313b1,
8.1771913681185391369911378989226b-22 %i
+ 1.2000000000000000000018070060519b1,
1.2999999999999999999970862549831b1 - 8.4503320198031190776637560809307b-22
%i, 7.1451883156443663949192698597447b-22 %i
+ 1.4000000000000000000035112730379b1,
1.499999999999999999996809992733b1 - 5.0695028726203639184930269555246b-22 %i,
1.6999999999999999999988916886292b1 - 1.4494744674447851086177248525705b-22
%i, 3.0226064149553728345607531928263b-22 %i
+ 1.6000000000000000000021831631434b1,
1.8999999999999999999999080171247b1 - 1.1346865485773691647618256857174b-23
%i, 5.0984652573716924251511605545068b-23 %i
+ 1.8000000000000000000003994624295b1,
1.1745417071012588556778686605423b-24 %i + 2.0000000000000000000000101739733b1]
Ray
/* -*- Mode: MACSYMA; Package: MAXIMA -*- */
/*
* Compute p/(x-r), returning p(r) and the new polynomial
* p is a list of the coefficients of the polynomial, arranged
* in descending powers.
*/
syndiv(p, r) :=
block([q : [], term : first(p)],
for c in rest(p) do
block([],
q : cons(term, q),
term : expand(rectform(r * term + c))
),
[term, reverse(q)]);
/* Evaluate polynomial at z */
poly_eval(a, z) :=
block([term : first(a)],
for c in rest(a) do
block([],
term : expand(rectform(z * term + c))
),
term);
/* Evaluate polynomial at z and also derivative of polynomial */
poly_eval_2(a, z) :=
block([q : 0, p : first(a)],
for c in rest(a) do
block([],
q : expand(rectform(p + z * q)),
p : expand(rectform(z * p + c))
),
[p, q]);
/* Evaluate polynomial at z, derivative of polynomial at z
* and an error bound. This is based on an algorithm given by Kahan
* http://www.cs.berkeley.edu/~wkahan/Math128/Poly.pdf.
*/
poly_eval_bound(a, z) :=
block([r : abs(z), q : 0, p : first(a), e : abs(first(a)) / 2],
for c in rest(a) do
block([],
q : expand(rectform(p + z * q)),
p : expand(rectform(z * p + c)),
e : expand(rectform(r * e + abs(p)))
),
e : expand(rectform((e - abs(p)) + e)),
[p, q, e]);
/*
* Compute a lower bound on the roots of the polynomial a
*
* If the polynomial a is sum a[k]*x^(n-k), k = 0,...,n, then this
* bound is the positive real root of the polynomial
*
* |a[0]|*x^n + |a[1]|*x^(n-1) + ... + |a[n-1]|*x - |a[n]|
*
* Use Newton-Raphson to find this root.
*/
lower_bound(a) :=
block([p : new_poly(a), root, f],
root : -p[length(p)]/p[length(p)-1],
do
(f : poly_eval_bound(p, root),
if abs(f[1]) < 2*f[3]*10^(-fpprec) then return (root),
root : root - f[1]/f[2]
)
);
/*
* Stage 1 shifts
*
* H(z) = [H(z) - H(0)/P(0)*P(z)]/z
*
* where we start with H(z) = P'(z), the derivative of the polynomial P.
*/
stage1(p, nshifts) :=
block([p0 : p[length(p)], deg : length(p) - 1, h, mult, maperror:false],
/* Compute derivative of p */
h : map(lambda([c,n], expand(c * (deg - n))),
p, makelist(i,i,0,deg-1)),
for k : 1 thru nshifts
do
(mult : expand(rectform(h[length(h)] / p0)),
/* 1/z*[H(z) - H(0)/P(0)*P(z)] */
h : map(lambda([hh,pp, nn], expand(hh - mult * pp)),
cons(0,h),
p,
makelist(i,i,1,length(p)-1))
),
h
);
/*
* Stage 2 shifts
*
* H(z) = [H(z) - H(s)/P(s)*P(z)]/(z-s)
*
* We start with H(z) computed from stage 1.
*
* s = b*exp(%i*phi), where b is the lower bound on the roots of P(z), and
* phi is some random phase.
*/
stage2(p, h, s, nshifts) :=
block([p0 : poly_eval(p, s),
t0 : 0,
t1 : 0,
t2, mult,
w, result : []],
for k : 1 thru nshifts do
(hv : poly_eval(h, s),
mult : expand(rectform(hv/p0)),
w : map(lambda([hh, pp], expand(hh - mult * pp)),
cons(0,h),
p),
t2 : expand(rectform(s - poly_eval(p, s) / (hv / h[1]))),
if (k >= 2) and (abs(t1 - t0) < abs(t0)/2) and (abs(t2 - t1) < abs(t1)/2) then
(result : [h, t2, k], return (result)),
t0 : t1, t1 : t2,
h : syndiv(w, s)[2]
),
result);
/* Determine if root is really a root of p.
*
* Return true if it is a root, and also the corresponding value of P(r).
*/
converged(p, root) :=
block([vals : poly_eval_bound(p, root)],
[is(abs(vals[1]) < 2*vals[3]*10^(-fpprec)), vals[1]]);
/*
* Stage 3 shifts
*
* Starting with the value of H(z) computed in stage 2, and
* s as the same value as in stage 2, perform the iteration:
*
* H(z) = [H(z) - H(s)/P(s)*P(z)]/(z-s)
*
* s = s - P(s)/(H(s)/h0)
*
* where h0 is the leading (highest power) coefficent of H(z).
*
* Perform this iteration until P(s) is as close to 0 as possible
* given roundoff errors.
*/
stage3(p, h, s, nshifts) :=
block([result : [], w, mult, root],
for k : 1 thru nshifts do
(root : s - expand(rectform(poly_eval(p,s)/(poly_eval(h,s)/h[1]))),
/*print("root = ", root),*/
vals : converged(p, root),
/*print("conv = ", vals),*/
if vals[1] then (result : [root, h, k], return (result)),
mult : expand(rectform(poly_eval(h, root) / vals[2])),
w : map(lambda([hh, pp], expand(hh - mult * pp)),
cons(0,h),
p),
h : syndiv(w, root)[2]),
result);
/*
* Main Jenkins-Traub iteration. Perform all three stages of iteration
*
* This needs more error checking in case one of the stages fails to converge.
* We should then either chose more iterations, or select a new shift s.
*/
jt_aux(p) :=
block([beta : lower_bound(p), s, h, h2, h3],
s : expand(rectform(beta * (cos(1.1) + %i*sin(1.1)))),
h : stage1(p, 5),
h2 : stage2(p, h, s, 50),
h3 : stage3(p, h2[1], s, 100),
h3[1]
);
/*
* Find one root of a polynomial P. It does the obvious thing if
* the degree of the polynomial is 2 or less. For higher degrees,
* Jenkins-Traub is used.
*/
find_one_root(p) :=
block([deg : length(p) - 1],
if is(degree = 0) then
[]
elseif is(degree = 1) then
-p[1]/p[2]
elseif is(degree = 2) then
block([a: p[1], b: p[2], c: p[3], s, discr],
s : if is(realpart(b) > 0) then 1 else -1,
discr : expand(rectform(sqrt(b*b - 4*a*c))),
(-b + s*discr)/(2*a)
)
else
jt_aux(p)
);
/*
* Find all roots of a polynomial p, if possible.
*
* P should be a list of the coefficients of the polynomial,
* arranged in descending powers. They should also be bfloats, but
* we don't check for that.
*
* A list of the roots, in the order in which they were computed.
* Multiple roots are listed multiple times. The roots are roughly in
* ascending order of magnitude.
*/
polyroots(p) :=
block([nroots : length(p)-1, roots : [], r],
for k : 1 thru nroots do
(r : find_one_root(p),
if r = [] then return(reverse(roots)),
roots: cons(r, roots),
p : syndiv(p, r)[2]),
reverse(roots));