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
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
00138 z = Complex(0);
00139 spir = Complex(-ONEPQT / 10, 1);
00140 do {
00141 pz = evalAndDerive(p, z, ppz, pppz);
00142
00143 zs = local::fejer(m, pz, ppz, pppz);
00144 f = abs(zs);
00145
00146
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
00164 w = sqm * ml0;
00165
00166
00167 ub = std::min(g, BIG_ONE * std::min(f, w));
00168 c = cauchyLowerBound(mP, ub);
00169
00170
00171 ub = std::min(RCONST * m * c, ub);
00172 lb = SML_ONE * c;
00173
00174
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
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
00200
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
00213
00214
00215
00216
00217 mpz0 = mpz;
00218 pz0 = pz;
00219 z0 = z;
00220 startd = 1;
00221
00222 }
00223 }
00224 }
00225 else {
00226
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
00256 zs = local::fejer(m, pz, ppz, pppz);
00257 f = abs(zs);
00258
00259
00260 step = local::laguerreStep(m, zs, pz, ppz);
00261 mstep = abs(step);
00262
00263
00264 w = sqm * mstep;
00265 f = std::min(w, f);
00266
00267
00268 ratio = mstep / g;
00269
00270 if (mstep < fpSpecs.ARE * abs(z))
00271 break;
00272 }
00273 }
00274
00275
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
00314 modulus(p, mP);
00315 mP[0] = -mP[0];
00316 }
00317 sortZeros(zeros);
00318 return (int)zeros.size();
00319 }
00320
00321 #endif //LAGUERRE_H