laguerre.h

Go to the documentation of this file.
00001 #ifndef LAGUERRE_H
00002 #define LAGUERRE_H
00003 
00029 #include "polyzero.h"
00030 
00061 template<class T, class U>
00062 int laguerreZeros(const Polynomial< T >& P, std::vector< std::complex < U > >& zeros, bool polish) {
00063 
00064     typedef U Float;
00065     typedef std::complex<Float> Complex;
00066     typedef Polynomial<Complex> ComplexPolynomial;
00067     typedef Polynomial<Float> FloatPolynomial;
00068 
00069     // local scope functions 
00070     class local {
00071     public:
00072         static Complex fejer(Float m, const Complex& pz, const Complex& ppz, const Complex& pppz) {
00073             Complex t1, t2, t3, z1, z2;
00074             t1 = pppz / (m * (m - 1));
00075             t2 = Float(2) * ppz / m;
00076             t3 = pz;
00077             solveDegree2(t1, t2, t3, z1, z2);
00078             return z1;
00079         }
00080 
00081         static Complex laguerreStep(Float m, const Complex& zs, const Complex& pz, const Complex& ppz) {
00082             Complex t;
00083             t = ((m - 2) * ppz) / (m * pz);
00084             t = (zs * t) + (m - 1);
00085             return zs / t;
00086         }
00087     };
00088 
00089     static FloatSpecs<Float> fpSpecs;
00090     const Float SMALL = Float(1.e-3);
00091     const Float BIG_ONE = Float(1.0001);
00092     const Float SML_ONE = Float(0.9999);
00093     const Float RCONST = Float(1.445);
00094     const Float ONEPQT = Float(1.25);
00095     const Float GAMA = Float(1);
00096     const Float THETA = Float(2);
00097 
00098     Polynomial<T> p = P;
00099     ComplexPolynomial q, q2;
00100     FloatPolynomial mP, d, p2, rem;
00101     int startd, iter, spiral, n;
00102     Complex spir, temp, r;
00103     Float sqm, m, ratio;
00104 
00105     Complex zs, l0, step, l, z0, z, pz, ppz, pppz, pz0;
00106     Float c, f, g, w, ml0, ub, lb, mstep, ml, mpz, mpz0;
00107 
00108     d.resize(3);
00109 
00110     modulus(p, mP);
00111     scalePoly(p, mP);
00112     mP[0] = -mP[0];
00113     
00114     n = removeNullZeros(p);
00115     zeros.assign(n, Complex(0));
00116 
00117     while (p.degree() > 0) {
00118 
00119         m = Float(p.degree());
00120 
00121         if (p.degree() == 1) {
00122             z = solveDegree1(p[1], p[0]);
00123             zeros.push_back(z);
00124             break;
00125         }
00126         if (p.degree() == 2) {
00127             solveDegree2(p[2], p[1], p[0], z0, z);
00128             if (polish) {
00129                 newtonZero(P, z, pz, mpz, false);
00130                 newtonZero(P, z0, pz, mpz, false);
00131             }
00132             zeros.push_back(z0);
00133             zeros.push_back(z);
00134             break;
00135         }
00136 
00137         //  get an initial estimate z, p(z)
00138         z = Complex(0);
00139         spir = Complex(-ONEPQT / 10, 1);
00140         do {
00141             pz = evalAndDerive(p, z, ppz, pppz);
00142             //  compute zs, fejer bound
00143             zs = local::fejer(m, pz, ppz, pppz);
00144             f = abs(zs);
00145 
00146             //  diro = l0
00147             l0 = local::laguerreStep(m, zs, p[0], p[1]);
00148             ml0 = abs(l0);
00149 
00150             if (_isnan(ml0)) {
00151                 if (z == Complex(0))
00152                     z = -p[0];
00153                 else
00154                     z = spir * z;
00155             }
00156             else
00157                 break;
00158         } while (1);
00159 
00160         sqm = sqrt(m);
00161         g = zerosGeometricMean(p);
00162 
00163         //  Laguerre bound w = sqrt(m).|step|;
00164         w = sqm * ml0;
00165 
00166         //  get cauchy bound
00167         ub = std::min(g, BIG_ONE * std::min(f, w));
00168         c = cauchyLowerBound(mP, ub);
00169 
00170         //  compute upper and lower bound of smallest zero
00171         ub = std::min(RCONST * m * c, ub);
00172         lb = SML_ONE * c;
00173 
00174         //  init first pass
00175         f = ub;
00176         g = ub;
00177         step = l0;
00178         mstep = ml0;
00179         ratio = mstep / g;
00180         mpz = abs(pz);
00181         mpz0 = mpz;
00182         spiral = 0;
00183         startd = 0;
00184         iter = 0;
00185         while (1) {
00186             ++iter;
00187 
00188             if (ratio > THETA) {
00189                 if (startd) {
00190                     // we've lost convergence, reduce step by half
00191                     ml *= Float(.5);
00192                     l *= Float(.5);
00193                     if (ml > fpSpecs.ARE * abs(z))
00194                         z = z0 + l;
00195                     else
00196                         break;
00197                 }
00198                 else {
00199                     // no convergence yet, try to catch a dip with an outward 
00200                     // spiral
00201                     if (spiral) {
00202                         z = spir * z;
00203                     }
00204                     else {
00205                         spiral = 1;
00206                         spir = Complex(-ONEPQT / m, 1);
00207                         ml = lb / (m * m);
00208                         z =  (step / mstep) * lb;
00209                     }
00210 
00211                     if (iter >= 25) {
00212                         //  we're spiraling too much, might as well let laguerre run its 
00213                         //  course (<40 iters), or we could be going for 100ks 
00214                         //  iterations !
00215                         //  CHECK: spiraling because the lower bound is too small ?
00216 
00217                         mpz0 = mpz;
00218                         pz0 = pz;
00219                         z0 = z;
00220                         startd = 1;
00221                         //printf("X");
00222                     }
00223                 }
00224             }
00225             else {
00226                 // converging, limit the step to keep converging
00227                 startd = 1;
00228                 if (ratio > GAMA && (spiral || lb <= g)) {
00229                     ratio = GAMA / ratio;
00230                     step /= ratio;
00231                     mstep /= ratio;
00232                 }
00233                 g = f;
00234                 l = step;
00235                 ml = mstep;
00236                 z0 = z;
00237                 mpz0 = mpz;
00238                 z += l;
00239             }
00240             _clearfp();
00241             
00242             pz = evalDeriveAndDeflate(p, z, ppz, pppz, q);
00243             
00244             mpz = abs(pz);
00245             if (_isnan(mpz))
00246                 mpz = fpSpecs.INFINY;
00247 
00248             if (mpz == 0)
00249                 break;
00250             
00251             if (mpz >= mpz0 && startd) {
00252                 ratio = THETA * BIG_ONE;
00253             }
00254             else {
00255                 // get Fejer bound
00256                 zs = local::fejer(m, pz, ppz, pppz);
00257                 f = abs(zs);
00258                 
00259                 //  compute next laguerre step
00260                 step = local::laguerreStep(m, zs, pz, ppz);
00261                 mstep = abs(step);
00262 
00263                 // Laguerre bound
00264                 w = sqm * mstep;
00265                 f = std::min(w, f);
00266 
00267                 //  g is Fejer bound from preceding point (z0)
00268                 ratio = mstep / g;
00269 
00270                 if (mstep < fpSpecs.ARE * abs(z))
00271                     break;
00272             }
00273         }
00274         // store the zero, different strategies for real and complex-valued
00275         // coeficients
00276         if (sizeof(T) == sizeof(U)) {
00277             if (abs(z.imag()) > abs(z.real()) * 20 * fpSpecs.MRE) {
00278                 if (polish) {
00279                     newtonZero(P, z, pz, mpz);
00280                     evalAndDeflate(p, z, q);
00281                 }
00282                 zeros.push_back(z);
00283                 z = conj(z);
00284                 if (polish)
00285                     newtonZero(P, z, pz, mpz);
00286                 evalAndDeflate(q, z, q2);
00287                 zeros.push_back(z);
00288                 n = q2.degree();
00289                 for (int i = 0; i <= n; ++i)
00290                     *(Float*)&(p[i]) = q2[i].real();
00291                 p.resize(n + 1);
00292             }
00293             else {
00294                 if (polish) {
00295                     newtonZero(P, z, pz, mpz);
00296                     evalAndDeflate(p, z, q);
00297                 }
00298                 zeros.push_back(z);
00299                 n = q.degree();
00300                 for (int i = 0; i <= n; ++i)
00301                     *(Float*)&(p[i]) = q[i].real();
00302                 p.resize(n + 1);
00303             }
00304         }
00305         else {
00306             if (polish) {
00307                 newtonZero(P, z, pz, mpz);
00308                 evalAndDeflate(p, z, q);
00309             }
00310             zeros.push_back(z);
00311             *(ComplexPolynomial*)&p = q;
00312         }
00313         //printf("%d\n", iter);
00314         modulus(p, mP);
00315         mP[0] = -mP[0];
00316     }
00317     sortZeros(zeros);
00318     return (int)zeros.size();
00319 }
00320 
00321 #endif //LAGUERRE_H

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