archive/share/trash/jtroot3.mac

   1 /* -*- Mode: MACSYMA; Package: MAXIMA -*- */
   2
   3 /*
   4  * This is an implementation of Jenkins-Traub polynomial root finder, as
   5  * described in "A Three-Stage Variable-Shift Iteration for Polynomial
   6  * Zeros and Its Relation to Generalized Rayleigh Iteraion", Numer. Math.,
   7  * 14, 252-263 (1970).
   8  *
   9  * http://dz-srv1.sub.uni-goettingen.de/sub/digbib/loader?did=D196932
  10  *
  11  * I, Raymond Toy, the author, hereby place this in the public domain.
  12  */
  13
  14 /* Max number of iterations to use in stage2 and stage3 iterations */
  15 polyroot_max_iterations : 20$
  16
  17 /* Max number of major passes */
  18 polyroot_max_passes : 2$
  19
  20 /*
  21  * Compute p/(x-r), returning p(r) and the new polynomial
  22  * p is a list of the coefficients of the polynomial, arranged
  23  * in descending powers.
  24  */
  25
  26 prt_syndiv(p, r) :=
  27   block([q : [], term : first(p)],
  28     for c in rest(p) do
  29       block([],
  30         q : cons(term, q),
  31         term : expand(rectform(r * term + c))
  32         ),
  33     [term, reverse(q)])$
  34
  35 /* Evaluate polynomial at z */
  36 prt_eval(a, z) :=
  37   block([term : first(a)],
  38     for c in rest(a) do
  39       block([],
  40         term : expand(rectform(z * term + c))
  41         ),
  42     term)$
  43
  44 /* Evaluate polynomial at z and also derivative of polynomial */
  45 prt_eval_2(a, z) :=
  46   block([q : 0, p : first(a)],
  47     for c in rest(a) do
  48       block([],
  49         q : expand(rectform(p + z * q)),
  50         p : expand(rectform(z * p + c))
  51         ),
  52     [p, q])$
  53
  54 /* Evaluate polynomial at z, derivative of polynomial at z
  55  * and an error bound.  This is based on an algorithm given by Kahan
  56  * http://www.cs.berkeley.edu/~wkahan/Math128/Poly.pdf.
  57  */
  58 prt_eval_bound(a, z) :=
  59   block([r : abs(z), q : 0, p : first(a), e : abs(first(a)) / 2],
  60     for c in rest(a) do
  61       block([],
  62         q : expand(rectform(p + z * q)),
  63         p : expand(rectform(z * p + c)),
  64         e : expand(rectform(r * e + abs(p)))
  65         ),
  66     e : expand(rectform((e - abs(p)) + e)),
  67     [p, q, e])$
  68
  69 /*
  70  * Given a polynomial P in X, extract the coefficients of the polynomial,
  71  * in descending powers of X.
  72 */
  73
  74 prt_coeff_list(p,x) :=
  75   if ratp(p) then
  76     block([keepfloat:true],
  77           p: ratnumer(rat(p,x)),
  78           reverse(makelist(prt_check_coeff(ratcoef(p,x,k)),k,0,hipow(p,x))))
  79   else
  80     (p : expand(p),
  81       reverse(makelist(prt_check_coeff(coeff(p,x,k)),k,0,hipow(p,x))))$
  82
  83 prt_check_coeff(c):=
  84   block([newc],
  85     newc: bfloat(c),
  86     if not(numberp(realpart(newc)) and numberp(imagpart(newc))) then
  87        error("poly_root can't handle non-numeric coefficients: ",c),
  88     newc)$
  89
  90 /* Compute estimates on the accuracy of the roots */
  91 prt_root_error1(a, z) :=
  92   block([r : abs(z), q : 0, p : first(a), e : abs(first(a)), d],
  93     if is(z - 0b0 = 0b0)
  94     then
  95       false
  96     else
  97       (d : -e/r,
  98         for c in rest(a) do
  99         /* Is expand actually needed below? */
 100         block([],
 101           q : expand(rectform(p + z * q)),
 102           d : expand(rectform(r*d + e + abs(q + q) - abs(p))),
 103           p : expand(rectform(z * p + c)),
 104           e : expand(rectform(r * e + abs(p)))
 105           ),
 106         e : expand(rectform(e - abs(p))),
 107         d : expand(rectform(d - abs(q))),
 108         [p, q, e, d]))$
 109
 110 /* Let roots be a list of the (estimated) roots of the polynomial
 111  * A in X.  Return a list of the estimated error in the roots.
 112  * If the error cannot be estimated for a particular root, set the error
 113  * to be false.
 114 */
 115 prt_root_error(x, a, roots) :=
 116   block([p : prt_coeff_list(a, x)],
 117     map(lambda([z],
 118         if is(z - 0b0 = 0b0)
 119         then
 120           false
 121         else
 122           block([bounds : prt_root_error1(p, z),
 123                  eps : 2*10b0^(-fpprec),
 124                  err],
 125              err : (abs(bounds[1]) + eps*bounds[3]) / (abs(bounds[2]) - eps*bounds[4]),
 126              if is(err > 0) then err else false)),
 127          roots))$
 128
 129 prt_new_poly(a) :=
 130   block([p : map(lambda([c], expand(rectform(abs(c)))), a)],
 131     p[length(p)] : -p[length(p)],
 132     p)$
 133
 134 /*
 135  * Compute a lower bound on the roots of the polynomial a
 136  *
 137  * If the polynomial a is sum a[k]*x^(n-k), k = 0,...,n, then this
 138  * bound is the positive real root of the polynomial
 139  *
 140  * |a[0]|*x^n + |a[1]|*x^(n-1) + ... + |a[n-1]|*x - |a[n]|
 141  *
 142  * Use Newton-Raphson to find this root.
 143  */
 144 prt_lower_bound(a) :=
 145   block([p : prt_new_poly(a), root, f],
 146     /*
 147      * Our initial guess is (|a[n]|/|a[0]|)^(1/n).  We should probably do
 148      * something better than that.
 149      */
 150     root : (-p[length(p)]/p[1])^(1/length(p)),
 151     do
 152       (f : prt_eval_bound(p, root),
 153         /* We don't need really good accuracy for the root here.*/
 154         if is(abs(f[1]) < 2*f[3]*10b0^(-fpprec/4)) then return (root),
 155         root : root - f[1]/f[2]
 156         )
 157       )$
 158
 159 /*
 160  * Stage 1 shifts
 161  *
 162  * H(z) = [H(z) - H(0)/P(0)*P(z)]/z
 163  *
 164  * where we start with H(z) = P'(z), the derivative of the polynomial P.
 165  */
 166
 167 prt_stage1(p, nshifts) :=
 168   block([p0 : p[length(p)], deg : length(p) - 1, h, mult,
 169          /* We want map to truncate silently, and without errors */
 170          mapprint : false, maperror : false],
 171     /* Compute derivative of p */
 172     h : map(lambda([c,n], expand(c * (deg - n))),
 173                     p, makelist(i,i,0,deg-1)),
 174     for k : 1 thru nshifts
 175       do
 176         (mult : expand(rectform(h[length(h)] / p0)),
 177          /* 1/z*[H(z) - H(0)/P(0)*P(z)] */
 178          h : map(lambda([hh,pp, nn], expand(hh - mult * pp)),
 179                  cons(0,h),
 180                  p,
 181                  makelist(i,i,1,length(p)-1))
 182          ),
 183        h
 184        )$
 185
 186 /*
 187  * Stage 2 shifts
 188  *
 189  * H(z) = [H(z) - H(s)/P(s)*P(z)]/(z-s)
 190  *
 191  * We start with H(z) computed from prt_stage 1.
 192  *
 193  */
 194 prt_stage2(p, h, s, nshifts) :=
 195   block([p0 : prt_eval(p, s),
 196          t0 : 0,
 197          t1 : 0,
 198          t2 : 0, mult, hv,
 199          w, result : []],
 200     for k : 1 thru nshifts do
 201       (hv : prt_eval(h, s),
 202         mult : expand(rectform(hv/p0)),
 203         w : map(lambda([hh, pp], expand(hh - mult * pp)),
 204                 cons(0,h),
 205                 p),
 206         /*
 207          * In some cases the leading coefficient of H can be zero, but this
 208          * only seems to happen on the first iteration, because stage1 can
 209          * produce such polynomials.  So skip it.
 210          */
 211         if (k > 1) then
 212           t2 : s - expand(rectform(prt_eval(p, s) / (hv / h[1]))),
 213         /*
 214          * t0, t1, t2 are successive estimates of the root.  We can
 215          * stop stage 2 when the differences between the roots are small enough.
 216          */
 217         if (k >= 2) and (abs(t1 - t0) < abs(t0)/2)
 218             and (abs(t2 - t1) < abs(t1)/2) then
 219           (result : [h, t2, k], return (result)),
 220         t0 : t1, t1 : t2,
 221         h : prt_syndiv(w, s)[2]
 222         ),
 223       result)$
 224
 225 /* Determine if root is really a root of p.
 226  *
 227  * Return true if it is a root, and also the corresponding value of P(r).
 228  */
 229 prt_converged(p, root) :=
 230   block([vals : prt_eval_bound(p, root)],
 231     [is(abs(vals[1]) < 2*vals[3]*(10*10b0^(-fpprec))), vals[1]])$
 232
 233 /*
 234  * Stage 3 shifts
 235  *
 236  * Starting with the value of H(z) computed in stage 2, and
 237  * s as the same value as in stage 2, perform the iteration:
 238  *
 239  * H(z) = [H(z) - H(s)/P(s)*P(z)]/(z-s)
 240  *
 241  * s = s - P(s)/(H(s)/h0)
 242  *
 243  * where h0 is the leading (highest power) coefficent of H(z).
 244  *
 245  * Perform this iteration until P(s) is as close to 0 as possible
 246  * given roundoff errors.
 247  */
 248 prt_stage3(p, h, s, nshifts) :=
 249   block([result : [], w, mult, root : s, vals],
 250     for k : 1 thru nshifts do
 251       (root : root - expand(rectform(prt_eval(p,root)/(prt_eval(h,root)/h[1]))),
 252        /*print("root = ", root),*/
 253        vals : prt_converged(p, root),
 254        /*print("conv = ", vals),*/
 255        if vals[1] then (result : [root, h, k], return (result)),
 256        mult : expand(rectform(prt_eval(h, root) / vals[2])),
 257        w : map(lambda([hh, pp], expand(hh - mult * pp)),
 258                cons(0,h),
 259                p),
 260        h : prt_syndiv(w, root)[2]),
 261      result)$
 262
 263 /*
 264  * Main Jenkins-Traub iteration.  Perform all three stages of iteration
 265  *
 266  * This needs more error checking in case one of the stages fails to converge.
 267  * We should then either chose more iterations, or select a new shift s.
 268  */
 269 prt_aux(p) :=
 270   block([beta : prt_lower_bound(p),
 271          /* r is a 94 deg rotation that we use for the desired shift */
 272          r : bfloat(rectform(exp(%i*%pi*94/180))),
 273          s, h, h2, h3, result : []],
 274     s : beta,
 275     /* Perform polyroot_max_passes major passes with different sequences */
 276     for pass : 1 thru polyroot_max_passes do
 277       for k : 1 thru 9 do
 278         (/* Initial shift is a 94 deg rotation of the previous shift */
 279          s : expand(rectform(s * r)),
 280          /* Perform 5 stage 1 shifts, then stage2 shifts */
 281          h : prt_stage1(p, 5),
 282          /*
 283           * I think the number of iterations for prt_stage2 and 3 might
 284           * depend on the desired precision.  Don't know how though.
 285           */
 286          h2 : prt_stage2(p, h, s, polyroot_max_iterations*pass),
 287          /* If stage 2 converged, try stage 3 */
 288          if h2 # [] then
 289            (h3 : prt_stage3(p, h2[1], s, polyroot_max_iterations),
 290             if h3 # [] then (result : h3[1], return(result)))),
 291   result)$
 292
 293 /*
 294  * Find one root of a polynomial P.  It does the obvious thing if
 295  * the degree of the polynomial is 2 or less.  For higher degrees,
 296  * Jenkins-Traub is used.
 297  */
 298 prt_find_one_root(p) :=
 299   block([degree : length(p) - 1],
 300     if is(degree = 0) then
 301       []
 302     elseif is(p[degree+1] = 0b0) then
 303       0
 304     elseif is(degree = 1) then
 305       rectform(-p[2]/p[1])
 306     elseif is(degree = 2) then
 307       block([a: p[1], b: p[2], c: p[3], s, discr],
 308         /*
 309          * Use the quadratic formula to find one root.  Try to minimize
 310          * roundoff by choosing the appropriate sign for the discriminant.
 311          */
 312         s : if is(realpart(b) > 0) then -1 else 1,
 313         discr : bfloat(expand(rectform(sqrt(b*b - 4*a*c)))),
 314         expand(rectform((-b + s*discr)/(2*a)))
 315         )
 316     else
 317       prt_aux(p)
 318       )$
 319
 320 /*
 321  * Find all roots (real and complex) of a polynomial p with
 322  * (real or complex) numerical coefficients.
 323  *
 324  * Returns a list of the roots, in the order in which they were computed.
 325  * Multiple roots are listed multiple times.  The roots are roughly in
 326  * ascending order of magnitude.
 327  *
 328  * The number of roots returned may not be the degree of the polynomial
 329  * if the algorithm failed to converge at any point.
 330  */
 331 polyroots(p, x) :=
 332   block([roots : [], nroots,
 333          r,
 334          poly : prt_coeff_list(p, x)],
 335     nroots : length(poly) - 1,
 336     for k : 1 thru nroots do
 337       (r : prt_find_one_root(poly),
 338        if r = [] then return(reverse(roots)),
 339        roots: cons(r, roots),
 340        poly : prt_syndiv(poly, r)[2]),
 341     reverse(roots))$
 342
 343 /* Here are some simple prt_tests, from TOMS Algorithm 419 */
 344 /* These should check their results for correctness! */
 345
 346 prt_test1() :=
 347   block([roots : sort(polyroots(product(x-k,k,1,10),x), "<")],
 348     print("Polynomial with roots 1,2,...,10"),
 349     roots)$
 350
 351 prt_test2() :=
 352   block([roots : sort(polyroots((x-%i)*(x-10000*%i)*(x-1/10000*%i),x), "<")],
 353     print("Polynomial with 3 roots on the imaginary axis"),
 354     roots)$
 355
 356 prt_test3() :=
 357   block([p : product(x-(1+%i)/2^k,k,0,9), roots],
 358     print("Polyomial with roots at (1+%i)/2^k, k = 0 to 9"),
 359     roots : sort(polyroots(p,x), "<"))$
 360
 361 prt_test4() :=
 362   block([p : (x-1)^4*(x-2*%i)^3*(x-3)^2*(x-4*%i), roots],
 363     print("Polynomial with repeated roots at 1, 2*%i, and 3"),
 364     roots : sort(polyroots(p,x), "<"))$
 365
 366 prt_test5() :=
 367   block([p : ((x-2*%i)^12+1), roots],
 368     print("Polynomial with roots evenly spaced on the unit circle centered at 2*%i"),
 369     roots : sort(polyroots(p,x), "<"))$