jenkinstraub.h

Go to the documentation of this file.
00001 #ifndef JENKINSTRAUB_H
00002 #define JENKINSTRAUB_H
00003 
00029 #include "polyzero.h"
00030 
00054 template<class T, class U>
00055 int jenkinsTraubZeros(const Polynomial< T >& P, std::vector< std::complex < U > >& zeros, bool polish) {
00056     
00057     typedef U Float;
00058     typedef std::complex<Float> Complex;
00059     typedef Polynomial<Complex> ComplexPolynomial;
00060     typedef Polynomial<Float> FloatPolynomial;
00061 
00062     class local {
00063     public:
00064         static void noShift(const ComplexPolynomial& p, ComplexPolynomial& h) {
00065             const unsigned M = 5;
00066             int m = p.degree();
00067             Float ii = Float(1);
00068             Float dd = Float(m);
00069             h.resize(m);
00070             //  first evalutate H0 = p'
00071             for (int i = 0; i < m; ++i, ii += Float(1))
00072                 h[i] = p[i + 1] * (ii / dd);
00073 
00074             // Hl+1(z) = 1/z * (Hl(z) - (Hl(0)/p(0)) * p(z));
00075             for (unsigned i = 0; i < M; ++i) {
00076                 Complex k = -h[0] / p[0];
00077                 for (int j = 0; j < m - 1; ++j)
00078                     h[j] = h[j + 1] + (k * p[j + 1]);
00079                 h[m - 1] = k * p[m];
00080             }
00081         }
00082 
00083         static bool fixedShift(unsigned pass, const ComplexPolynomial& p, ComplexPolynomial& h, ComplexPolynomial& qp, Complex& s) {
00084 
00085             bool failed = false;
00086             int conv = 0;
00087             Complex hs, ps, t, tOld, z;
00088             ComplexPolynomial qh;
00089             ps = p.evalAndDeflate(s, qp);
00090             t = nextT(s, ps, h, hs, qh);
00091 
00092             for (unsigned i = 0; i < (10 * pass); ++i) {
00093                 tOld = t;
00094                 //  compute next h and new t
00095                 nextH(p, h, qp, qh, hs, ps, t);
00096                 t = nextT(s, ps, h, hs, qh);
00097                 z = s + t;
00098 
00099                 if (hs == Complex(0))
00100                     continue;
00101 
00102                 // weak convergence test
00103                 Float k = abs(tOld - t);
00104                 Float l = Float(.5) * abs(z);
00105                 if (!failed && abs(t - tOld) < Float(.5) * abs(z)) {
00106                     if (++conv >= 2) {
00107                         if (variableShift(p, h, qp, qh, z, t, hs, ps)) {
00108                             s = z;
00109                             return true;
00110                         }
00111                         conv = 0;
00112                         failed = true;
00113                     }
00114                 }
00115                 else
00116                     conv = 0;
00117             }
00118             // do a last attempt with final h polynomial
00119             if (variableShift(p, h, qp, qh, z, t, ps, hs)) {
00120                 s = z;
00121                 return true;
00122             }
00123             return false;
00124         }
00125 
00126         static bool variableShift(const ComplexPolynomial& p, ComplexPolynomial& h, ComplexPolynomial& qp, ComplexPolynomial& qh, Complex& s, Complex& t, Complex& ps, Complex& hs) {
00127             static FloatSpecs<Float> fpSpecs;
00128             bool b = false;
00129             int i;
00130             Float step = 0, mpOld = fpSpecs.INFINY, mp = 0, mt = fpSpecs.INFINY, mtOld = fpSpecs.INFINY;
00131             ComplexPolynomial _H = h, _QH = qh;
00132             Complex _s = s, _Hs = hs, _t = t;
00133             Float e;
00134 
00135             for (i = 1; i <= 12; ++i) {
00136                 ps = evalAndDeflate(p, s, qp, e);
00137                 mpOld = mp;
00138                 mp = abs(ps);
00139                 mtOld = mt;
00140                 mt = abs(t);
00141 
00142                 if (_isnan(mp)) {
00143                     //printf("a");
00144                     break;
00145                 }
00146 
00147                 if (mp < e) {
00148                     //printf("i:%d ", i);
00149                     return true;
00150                 }
00151 
00152                 if (i != 0) {
00153                     if(!b && mp > mpOld && step < Float(.05)) {
00154                         // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed 
00155                         // shift steps into the cluster to force one zero to dominate
00156                         Float tp = std::max(step, fpSpecs.ETA);
00157                         b = true;
00158                         Complex r(1 + sqrt(tp), sqrt(tp));
00159                         s *= r;
00160                         ps = p.evalAndDeflate(s, qp);
00161                         for(unsigned j = 0; j < 5; ++j) {
00162                             t = nextT(s, ps, h, hs, qh);
00163                             nextH(p, h, qp, qh, hs, ps, t);
00164                         }
00165                         mp = fpSpecs.INFINY;
00166                     }
00167                     // Exit if polynomial value increase significantly
00168                     else if((mp * Float(0.1)) > mpOld) 
00169                         break;
00170                 }
00171 
00172                 t = nextT(s, ps, h, hs, qh);
00173                 nextH(p, h, qp, qh, hs, ps, t);
00174                 t = nextT(s, ps, h, hs, qh);
00175                 if (hs != Complex(0)) {
00176                     step = abs(t) / abs(s);
00177                     s = s + t;
00178                 }
00179             }
00180             h = _H;
00181             qh = _QH;
00182             s = _s;
00183             t = _t;
00184             hs = _Hs;
00185             ps = evalAndDeflate(p, s, qp);
00186             //printf("i:%d ", i);
00187             return false;
00188         }
00189 
00190         static Complex nextT(const Complex& s, const Complex& ps, const ComplexPolynomial& h, Complex& hs, ComplexPolynomial& qh) {
00191             static FloatSpecs<Float> fpSpecs;
00192             Complex r;
00193             Float e;
00194             hs = evalAndDeflate(h, s, qh, e);
00195             if (abs(hs) < fpSpecs.ARE * (10 * abs(h[0]))) {
00196             //if (abs(hs) < e) {
00197                 hs = Complex(0);
00198                 return Complex(0);
00199             }
00200             return -ps / hs;
00201         }
00202 
00203         static void nextH(const ComplexPolynomial& p, ComplexPolynomial& h, const ComplexPolynomial& qp, const ComplexPolynomial& qh, const Complex& hs, const Complex& ps, const Complex& t) {
00204             int m = qp.degree();
00205             int n = qh.degree();
00206             int i;
00207             if (hs != Complex(0)) {
00208                 // Hl+1(z) = 1/(z-s) * Hl(z) - (Hl(s)/p(s)) * p(z))
00209                 h.resize(m + 1);
00210                 for(i = 0; i <= n; ++i)
00211                     h[i] = (t * qh[i]) + qp[i];
00212                 for( ; i <= m; ++i)
00213                     h[i] = qp[i];
00214             }
00215             else
00216                 h = qh;
00217         }
00218     };
00219 
00220     const bool ADAPTIVE = false;
00221     const int MAXIT1 = 2;
00222     const int MAXIT2 = 9;
00223     static const Complex SHIFT = exp(Complex(0, Float(-94) * acos(Float(-1)) / Float(180)));
00224 
00225     ComplexPolynomial p, qp, qp2, h;
00226     FloatPolynomial mP;
00227     Complex z, pz, z1, pz1;
00228     Float e, mpz;
00229     int i, n;
00230 
00231     p.resize(P.degree() + 1);
00232     for (i = 0; i < (int)p.size(); ++i)
00233         p[i] = Complex(P[i]);
00234 
00235     zeros.clear();
00236 
00237     modulus(p, mP);
00238     scalePoly(p, mP);
00239 
00240     i = removeNullZeros(p);
00241     zeros.assign(i, Complex(0));
00242 
00243     Complex shift = SHIFT;
00244     while (p.degree() > 0) {
00245 
00246         if (p.degree() == 1) {
00247             z = solveDegree1(p[1], p[0]);
00248             zeros.push_back(z);
00249             break;
00250         }
00251         if (p.degree() == 2) {
00252             solveDegree2(p[2], p[1], p[0], z, z1);
00253             zeros.push_back(z);
00254             zeros.push_back(z1);
00255             break;
00256         }
00257 
00258         modulus(p, mP);
00259         mP[0] = -mP[0];
00260         Float beta = cauchyLowerBound(mP) * Float(1.05);
00261         bool zeroFound = false;
00262 
00263         for (i = 1; i < MAXIT1 && !zeroFound; ++i) {
00264             // stage1, no shift
00265             local::noShift(p, h);
00266             for (unsigned j = 1; j < MAXIT2 && !zeroFound; ++j) {
00267                 // stage2, fixed shift
00268                 z = beta * shift;
00269                 shift *= SHIFT;
00270                 if (local::fixedShift(j, p, h, qp, z)) {
00271                     zeroFound = true;
00272                     break;
00273                 }
00274             }
00275         }
00276         if (zeroFound) {
00277             //printf("\n");
00278             if (sizeof(T) == sizeof(U)) {
00279                 z1 = conj(z);
00280                 pz1 = eval(qp, conj(z), e);
00281                 if (abs(pz1) < e) {
00282                     if (polish) {
00283                         newtonZero(P, z, pz, mpz, ADAPTIVE);
00284                         evalAndDeflate(p, z, qp);
00285                     }
00286                     zeros.push_back(z);
00287                     z = conj(z);
00288                     if (polish) {
00289                         newtonZero(P, z, pz, mpz, ADAPTIVE);
00290                     }
00291                     evalAndDeflate(qp, z, qp2);
00292                     zeros.push_back(z);
00293                     n = qp2.degree();
00294                     for (int i = 0; i <= n; ++i)
00295                         *(Float*)&(p[i]) = qp2[i].real();
00296                     p.resize(n + 1);
00297                 }
00298                 else {
00299                     if (polish) {
00300                         newtonZero(P, z, pz, mpz, ADAPTIVE);
00301                         evalAndDeflate(p, z, qp);
00302                     }
00303                     zeros.push_back(z);
00304                     n = qp.degree();
00305                     for (int i = 0; i <= n; ++i)
00306                         *(Float*)&(p[i]) = qp[i].real();
00307                     p.resize(n + 1);
00308                 }
00309             }
00310             else {
00311                 if (polish) {
00312                     newtonZero(P, z, pz, mpz, ADAPTIVE);
00313                     evalAndDeflate(p, z, qp);
00314                 }
00315                 zeros.push_back(z);
00316                 *(ComplexPolynomial*)&p = qp;
00317             }
00318         }
00319         else 
00320             break;
00321     }
00322     sortZeros(zeros);
00323     return (int)zeros.size();
00324 }
00325 
00326 #endif //JENKINSTRAUB_H

Generated on Mon Aug 21 21:00:37 2006 for The Polynomials Templates Library by  doxygen 1.4.5