a00105.html

This library defines basic operation on polynomials, and contains also 3 different roots (zeroes)-fi

<a name="l00183"></a>00183         startd = 0;
<a name="l00184"></a>00184         iter = 0;
<a name="l00185"></a>00185         <span class="keywordflow">while</span> (1) {
<a name="l00186"></a>00186             ++iter;
<a name="l00187"></a>00187 
<a name="l00188"></a>00188             <span class="keywordflow">if</span> (ratio &gt; THETA) {
<a name="l00189"></a>00189                 <span class="keywordflow">if</span> (startd) {
<a name="l00190"></a>00190                     <span class="comment">// we've lost convergence, reduce step by half</span>
<a name="l00191"></a>00191                     ml *= Float(.5);
<a name="l00192"></a>00192                     l *= Float(.5);
<a name="l00193"></a>00193                     <span class="keywordflow">if</span> (ml &gt; fpSpecs.ARE * abs(z))
<a name="l00194"></a>00194                         z = z0 + l;
<a name="l00195"></a>00195                     <span class="keywordflow">else</span>
<a name="l00196"></a>00196                         <span class="keywordflow">break</span>;
<a name="l00197"></a>00197                 }
<a name="l00198"></a>00198                 <span class="keywordflow">else</span> {
<a name="l00199"></a>00199                     <span class="comment">// no convergence yet, try to catch a dip with an outward </span>
<a name="l00200"></a>00200                     <span class="comment">// spiral</span>
<a name="l00201"></a>00201                     <span class="keywordflow">if</span> (spiral) {
<a name="l00202"></a>00202                         z = spir * z;
<a name="l00203"></a>00203                     }
<a name="l00204"></a>00204                     <span class="keywordflow">else</span> {
<a name="l00205"></a>00205                         spiral = 1;
<a name="l00206"></a>00206                         spir = Complex(-ONEPQT / m, 1);
<a name="l00207"></a>00207                         ml = lb / (m * m);
<a name="l00208"></a>00208                         z =  (step / mstep) * lb;
<a name="l00209"></a>00209                     }
<a name="l00210"></a>00210 
<a name="l00211"></a>00211                     <span class="keywordflow">if</span> (iter &gt;= 25) {
<a name="l00212"></a>00212                         <span class="comment">//  we're spiraling too much, might as well let laguerre run its </span>
<a name="l00213"></a>00213                         <span class="comment">//  course (&lt;40 iters), or we could be going for 100ks </span>
<a name="l00214"></a>00214                         <span class="comment">//  iterations !</span>
<a name="l00215"></a>00215                         <span class="comment">//  CHECK: spiraling because the lower bound is too small ?</span>
<a name="l00216"></a>00216 
<a name="l00217"></a>00217                         mpz0 = mpz;
<a name="l00218"></a>00218                         pz0 = pz;
<a name="l00219"></a>00219                         z0 = z;
<a name="l00220"></a>00220                         startd = 1;
<a name="l00221"></a>00221                         <span class="comment">//printf("X");</span>
<a name="l00222"></a>00222                     }
<a name="l00223"></a>00223                 }
<a name="l00224"></a>00224             }
<a name="l00225"></a>00225             <span class="keywordflow">else</span> {
<a name="l00226"></a>00226                 <span class="comment">// converging, limit the step to keep converging</span>
<a name="l00227"></a>00227                 startd = 1;
<a name="l00228"></a>00228                 <span class="keywordflow">if</span> (ratio &gt; GAMA &amp;&amp; (spiral || lb &lt;= g)) {
<a name="l00229"></a>00229                     ratio = GAMA / ratio;
<a name="l00230"></a>00230                     step /= ratio;
<a name="l00231"></a>00231                     mstep /= ratio;
<a name="l00232"></a>00232                 }
<a name="l00233"></a>00233                 g = f;
<a name="l00234"></a>00234                 l = step;
<a name="l00235"></a>00235                 ml = mstep;
<a name="l00236"></a>00236                 z0 = z;
<a name="l00237"></a>00237                 mpz0 = mpz;
<a name="l00238"></a>00238                 z += l;
<a name="l00239"></a>00239             }
<a name="l00240"></a>00240             _clearfp();
<a name="l00241"></a>00241             
<a name="l00242"></a>00242             pz = <a class="code" href="a00100.html#g055ddecc1576e178f623b4b8e868cc23">evalDeriveAndDeflate</a>(p, z, ppz, pppz, q);
<a name="l00243"></a>00243             
<a name="l00244"></a>00244             mpz = abs(pz);
<a name="l00245"></a>00245             <span class="keywordflow">if</span> (_isnan(mpz))
<a name="l00246"></a>00246                 mpz = fpSpecs.INFINY;
<a name="l00247"></a>00247 
<a name="l00248"></a>00248             <span class="keywordflow">if</span> (mpz == 0)
<a name="l00249"></a>00249                 <span class="keywordflow">break</span>;
<a name="l00250"></a>00250             
<a name="l00251"></a>00251             <span class="keywordflow">if</span> (mpz &gt;= mpz0 &amp;&amp; startd) {
<a name="l00252"></a>00252                 ratio = THETA * BIG_ONE;
<a name="l00253"></a>00253             }
<a name="l00254"></a>00254             <span class="keywordflow">else</span> {
<a name="l00255"></a>00255                 <span class="comment">// get Fejer bound</span>
<a name="l00256"></a>00256                 zs = local::fejer(m, pz, ppz, pppz);
<a name="l00257"></a>00257                 f = abs(zs);
<a name="l00258"></a>00258                 
<a name="l00259"></a>00259                 <span class="comment">//  compute next laguerre step</span>
<a name="l00260"></a>00260                 step = local::laguerreStep(m, zs, pz, ppz);
<a name="l00261"></a>00261                 mstep = abs(step);
<a name="l00262"></a>00262 
<a name="l00263"></a>00263                 <span class="comment">// Laguerre bound</span>
<a name="l00264"></a>00264                 w = sqm * mstep;
<a name="l00265"></a>00265                 f = std::min(w, f);
<a name="l00266"></a>00266 
<a name="l00267"></a>00267                 <span class="comment">//  g is Fejer bound from preceding point (z0)</span>
<a name="l00268"></a>00268                 ratio = mstep / g;
<a name="l00269"></a>00269 
<a name="l00270"></a>00270                 <span class="keywordflow">if</span> (mstep &lt; fpSpecs.ARE * abs(z))
<a name="l00271"></a>00271                     <span class="keywordflow">break</span>;
<a name="l00272"></a>00272             }
<a name="l00273"></a>00273         }
<a name="l00274"></a>00274         <span class="comment">// store the zero, different strategies for real and complex-valued</span>
<a name="l00275"></a>00275         <span class="comment">// coeficients</span>
<a name="l00276"></a>00276         <span class="keywordflow">if</span> (<span class="keyword">sizeof</span>(T) == <span class="keyword">sizeof</span>(U)) {
<a name="l00277"></a>00277             <span class="keywordflow">if</span> (abs(z.imag()) &gt; abs(z.real()) * 20 * fpSpecs.MRE) {
<a name="l00278"></a>00278                 <span class="keywordflow">if</span> (polish) {
<a name="l00279"></a>00279                     <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz);
<a name="l00280"></a>00280                     <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, q);
<a name="l00281"></a>00281                 }
<a name="l00282"></a>00282                 zeros.push_back(z);
<a name="l00283"></a>00283                 z = conj(z);
<a name="l00284"></a>00284                 <span class="keywordflow">if</span> (polish)
<a name="l00285"></a>00285                     <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz);
<a name="l00286"></a>00286                 <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(q, z, q2);
<a name="l00287"></a>00287                 zeros.push_back(z);
<a name="l00288"></a>00288                 n = q2.degree();
<a name="l00289"></a>00289                 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt;= n; ++i)
<a name="l00290"></a>00290                     *(Float*)&amp;(p[i]) = q2[i].real();
<a name="l00291"></a>00291                 p.resize(n + 1);
<a name="l00292"></a>00292             }
<a name="l00293"></a>00293             <span class="keywordflow">else</span> {
<a name="l00294"></a>00294                 <span class="keywordflow">if</span> (polish) {
<a name="l00295"></a>00295                     <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz);
<a name="l00296"></a>00296                     <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, q);
<a name="l00297"></a>00297                 }
<a name="l00298"></a>00298                 zeros.push_back(z);
<a name="l00299"></a>00299                 n = q.degree();
<a name="l00300"></a>00300                 <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt;= n; ++i)
<a name="l00301"></a>00301                     *(Float*)&amp;(p[i]) = q[i].real();
<a name="l00302"></a>00302                 p.resize(n + 1);
<a name="l00303"></a>00303             }
<a name="l00304"></a>00304         }
<a name="l00305"></a>00305         <span class="keywordflow">else</span> {
<a name="l00306"></a>00306             <span class="keywordflow">if</span> (polish) {
<a name="l00307"></a>00307                 <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz);
<a name="l00308"></a>00308                 <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, q);
<a name="l00309"></a>00309             }
<a name="l00310"></a>00310             zeros.push_back(z);
<a name="l00311"></a>00311             *(ComplexPolynomial*)&amp;p = q;
<a name="l00312"></a>00312         }
<a name="l00313"></a>00313         <span class="comment">//printf("%d\n", iter);</span>
<a name="l00314"></a>00314         <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mP);
<a name="l00315"></a>00315         mP[0] = -mP[0];
<a name="l00316"></a>00316     }
<a name="l00317"></a>00317     <a class="code" href="a00102.html#g3536daa88b2578a96cd84d2680ac9d66">sortZeros</a>(zeros);
<a name="l00318"></a>00318     <span class="keywordflow">return</span> (<span class="keywordtype">int</span>)zeros.size();
<a name="l00319"></a>00319 }
<a name="l00320"></a>00320 
<a name="l00321"></a>00321 <span class="preprocessor">#endif //LAGUERRE_H</span>
</pre></div><hr size="1"><address style="align: right;"><small>Generated on Mon Aug 21 21:57:25 2006 for The Polynomials Templates Library by&nbsp;
<a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border="0"></a> 1.4.5 </small></address>
</body>
</html>