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
00071 for (int i = 0; i < m; ++i, ii += Float(1))
00072 h[i] = p[i + 1] * (ii / dd);
00073
00074
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
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
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
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
00144 break;
00145 }
00146
00147 if (mp < e) {
00148
00149 return true;
00150 }
00151
00152 if (i != 0) {
00153 if(!b && mp > mpOld && step < Float(.05)) {
00154
00155
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
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
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
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
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
00265 local::noShift(p, h);
00266 for (unsigned j = 1; j < MAXIT2 && !zeroFound; ++j) {
00267
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
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