Hi
I'm trying to use your program for polynomial given by recurence
relation ( see: http://fraktal.republika.pl/mset_centers.html);
My polynomial :
P(n):=if n=0 then 0 else P(n-1)^2+c;
List of coefficients:
give_coefficients(n):=block( P(n):=if n=0 then 0 else P(n-1)^2+c,
Pn:expand(P(n)), degree:hipow(Pn,c),
a:makelist(coeff(Pn,c,degree-i),i,0,degree) );
example of use:
a:give_coefficients(3);
(%o7) [1,2,1,1,0]
but :
polyroots(a);
does not work.
What I do wrong?
Adam Majewski
Raymond Toy (RT/EUS) pisze:
> 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));