a00104.html

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

<a name="l00182"></a>00182             s = _s;
<a name="l00183"></a>00183             t = _t;
<a name="l00184"></a>00184             hs = _Hs;
<a name="l00185"></a>00185             ps = <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, s, qp);
<a name="l00186"></a>00186             <span class="comment">//printf("i:%d ", i);</span>
<a name="l00187"></a>00187             <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00188"></a>00188         }
<a name="l00189"></a>00189 
<a name="l00190"></a>00190         <span class="keyword">static</span> Complex nextT(<span class="keyword">const</span> Complex&amp; s, <span class="keyword">const</span> Complex&amp; ps, <span class="keyword">const</span> ComplexPolynomial&amp; h, Complex&amp; hs, ComplexPolynomial&amp; qh) {
<a name="l00191"></a>00191             <span class="keyword">static</span> FloatSpecs&lt;Float&gt; fpSpecs;
<a name="l00192"></a>00192             Complex r;
<a name="l00193"></a>00193             Float e;
<a name="l00194"></a>00194             hs = <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(h, s, qh, e);
<a name="l00195"></a>00195             <span class="keywordflow">if</span> (abs(hs) &lt; fpSpecs.<a class="code" href="a00089.html#6cf389757399a3d8def3269f8d1f7e81">ARE</a> * (10 * abs(h[0]))) {
<a name="l00196"></a>00196             <span class="comment">//if (abs(hs) &lt; e) {</span>
<a name="l00197"></a>00197                 hs = Complex(0);
<a name="l00198"></a>00198                 <span class="keywordflow">return</span> Complex(0);
<a name="l00199"></a>00199             }
<a name="l00200"></a>00200             <span class="keywordflow">return</span> -ps / hs;
<a name="l00201"></a>00201         }
<a name="l00202"></a>00202 
<a name="l00203"></a>00203         <span class="keyword">static</span> <span class="keywordtype">void</span> nextH(<span class="keyword">const</span> ComplexPolynomial&amp; p, ComplexPolynomial&amp; h, <span class="keyword">const</span> ComplexPolynomial&amp; qp, <span class="keyword">const</span> ComplexPolynomial&amp; qh, <span class="keyword">const</span> Complex&amp; hs, <span class="keyword">const</span> Complex&amp; ps, <span class="keyword">const</span> Complex&amp; t) {
<a name="l00204"></a>00204             <span class="keywordtype">int</span> m = qp.degree();
<a name="l00205"></a>00205             <span class="keywordtype">int</span> n = qh.degree();
<a name="l00206"></a>00206             <span class="keywordtype">int</span> i;
<a name="l00207"></a>00207             <span class="keywordflow">if</span> (hs != Complex(0)) {
<a name="l00208"></a>00208                 <span class="comment">// Hl+1(z) = 1/(z-s) * Hl(z) - (Hl(s)/p(s)) * p(z))</span>
<a name="l00209"></a>00209                 h.resize(m + 1);
<a name="l00210"></a>00210                 <span class="keywordflow">for</span>(i = 0; i &lt;= n; ++i)
<a name="l00211"></a>00211                     h[i] = (t * qh[i]) + qp[i];
<a name="l00212"></a>00212                 <span class="keywordflow">for</span>( ; i &lt;= m; ++i)
<a name="l00213"></a>00213                     h[i] = qp[i];
<a name="l00214"></a>00214             }
<a name="l00215"></a>00215             <span class="keywordflow">else</span>
<a name="l00216"></a>00216                 h = qh;
<a name="l00217"></a>00217         }
<a name="l00218"></a>00218     };
<a name="l00219"></a>00219 
<a name="l00220"></a>00220     <span class="keyword">const</span> <span class="keywordtype">bool</span> ADAPTIVE = <span class="keyword">false</span>;
<a name="l00221"></a>00221     <span class="keyword">const</span> <span class="keywordtype">int</span> MAXIT1 = 2;
<a name="l00222"></a>00222     <span class="keyword">const</span> <span class="keywordtype">int</span> MAXIT2 = 9;
<a name="l00223"></a>00223     <span class="keyword">static</span> <span class="keyword">const</span> Complex SHIFT = exp(Complex(0, Float(-94) * acos(Float(-1)) / Float(180)));
<a name="l00224"></a>00224 
<a name="l00225"></a>00225     ComplexPolynomial p, qp, qp2, h;
<a name="l00226"></a>00226     FloatPolynomial mP;
<a name="l00227"></a>00227     Complex z, pz, z1, pz1;
<a name="l00228"></a>00228     Float e, mpz;
<a name="l00229"></a>00229     <span class="keywordtype">int</span> i, n;
<a name="l00230"></a>00230 
<a name="l00231"></a>00231     p.resize(P.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>() + 1);
<a name="l00232"></a>00232     <span class="keywordflow">for</span> (i = 0; i &lt; (<span class="keywordtype">int</span>)p.size(); ++i)
<a name="l00233"></a>00233         p[i] = Complex(P[i]);
<a name="l00234"></a>00234 
<a name="l00235"></a>00235     zeros.clear();
<a name="l00236"></a>00236 
<a name="l00237"></a>00237     <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mP);
<a name="l00238"></a>00238     <a class="code" href="a00099.html#g663581aef6853ed86df03a35505c9ebd">scalePoly</a>(p, mP);
<a name="l00239"></a>00239 
<a name="l00240"></a>00240     i = <a class="code" href="a00102.html#g206394af98d861f1a81bea25ae6d84da">removeNullZeros</a>(p);
<a name="l00241"></a>00241     zeros.assign(i, Complex(0));
<a name="l00242"></a>00242 
<a name="l00243"></a>00243     Complex shift = SHIFT;
<a name="l00244"></a>00244     <span class="keywordflow">while</span> (p.degree() &gt; 0) {
<a name="l00245"></a>00245 
<a name="l00246"></a>00246         <span class="keywordflow">if</span> (p.degree() == 1) {
<a name="l00247"></a>00247             z = <a class="code" href="a00102.html#g6e9fe1cec8c7ec334d9e21f31b602bde">solveDegree1</a>(p[1], p[0]);
<a name="l00248"></a>00248             zeros.push_back(z);
<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> (p.degree() == 2) {
<a name="l00252"></a>00252             <a class="code" href="a00102.html#gb88e92200690d2e3f8b0c785c8fa42d3">solveDegree2</a>(p[2], p[1], p[0], z, z1);
<a name="l00253"></a>00253             zeros.push_back(z);
<a name="l00254"></a>00254             zeros.push_back(z1);
<a name="l00255"></a>00255             <span class="keywordflow">break</span>;
<a name="l00256"></a>00256         }
<a name="l00257"></a>00257 
<a name="l00258"></a>00258         <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mP);
<a name="l00259"></a>00259         mP[0] = -mP[0];
<a name="l00260"></a>00260         Float beta = <a class="code" href="a00099.html#g4c00b88d043d1f70cdcb4d2cec66b55e">cauchyLowerBound</a>(mP) * Float(1.05);
<a name="l00261"></a>00261         <span class="keywordtype">bool</span> zeroFound = <span class="keyword">false</span>;
<a name="l00262"></a>00262 
<a name="l00263"></a>00263         <span class="keywordflow">for</span> (i = 1; i &lt; MAXIT1 &amp;&amp; !zeroFound; ++i) {
<a name="l00264"></a>00264             <span class="comment">// stage1, no shift</span>
<a name="l00265"></a>00265             local::noShift(p, h);
<a name="l00266"></a>00266             <span class="keywordflow">for</span> (<span class="keywordtype">unsigned</span> j = 1; j &lt; MAXIT2 &amp;&amp; !zeroFound; ++j) {
<a name="l00267"></a>00267                 <span class="comment">// stage2, fixed shift</span>
<a name="l00268"></a>00268                 z = beta * shift;
<a name="l00269"></a>00269                 shift *= SHIFT;
<a name="l00270"></a>00270                 <span class="keywordflow">if</span> (local::fixedShift(j, p, h, qp, z)) {
<a name="l00271"></a>00271                     zeroFound = <span class="keyword">true</span>;
<a name="l00272"></a>00272                     <span class="keywordflow">break</span>;
<a name="l00273"></a>00273                 }
<a name="l00274"></a>00274             }
<a name="l00275"></a>00275         }
<a name="l00276"></a>00276         <span class="keywordflow">if</span> (zeroFound) {
<a name="l00277"></a>00277             <span class="comment">//printf("\n");</span>
<a name="l00278"></a>00278             <span class="keywordflow">if</span> (<span class="keyword">sizeof</span>(T) == <span class="keyword">sizeof</span>(U)) {
<a name="l00279"></a>00279                 z1 = conj(z);
<a name="l00280"></a>00280                 pz1 = <a class="code" href="a00100.html#gc5161c8ecd75b0108d0ea60c92d30d21">eval</a>(qp, conj(z), e);
<a name="l00281"></a>00281                 <span class="keywordflow">if</span> (abs(pz1) &lt; e) {
<a name="l00282"></a>00282                     <span class="keywordflow">if</span> (polish) {
<a name="l00283"></a>00283                         <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, ADAPTIVE);
<a name="l00284"></a>00284                         <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, qp);
<a name="l00285"></a>00285                     }
<a name="l00286"></a>00286                     zeros.push_back(z);
<a name="l00287"></a>00287                     z = conj(z);
<a name="l00288"></a>00288                     <span class="keywordflow">if</span> (polish) {
<a name="l00289"></a>00289                         <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, ADAPTIVE);
<a name="l00290"></a>00290                     }
<a name="l00291"></a>00291                     <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(qp, z, qp2);
<a name="l00292"></a>00292                     zeros.push_back(z);
<a name="l00293"></a>00293                     n = qp2.degree();
<a name="l00294"></a>00294                     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt;= n; ++i)
<a name="l00295"></a>00295                         *(Float*)&amp;(p[i]) = qp2[i].real();
<a name="l00296"></a>00296                     p.resize(n + 1);
<a name="l00297"></a>00297                 }
<a name="l00298"></a>00298                 <span class="keywordflow">else</span> {
<a name="l00299"></a>00299                     <span class="keywordflow">if</span> (polish) {
<a name="l00300"></a>00300                         <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, ADAPTIVE);
<a name="l00301"></a>00301                         <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, qp);
<a name="l00302"></a>00302                     }
<a name="l00303"></a>00303                     zeros.push_back(z);
<a name="l00304"></a>00304                     n = qp.degree();
<a name="l00305"></a>00305                     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt;= n; ++i)
<a name="l00306"></a>00306                         *(Float*)&amp;(p[i]) = qp[i].real();
<a name="l00307"></a>00307                     p.resize(n + 1);
<a name="l00308"></a>00308                 }
<a name="l00309"></a>00309             }
<a name="l00310"></a>00310             <span class="keywordflow">else</span> {
<a name="l00311"></a>00311                 <span class="keywordflow">if</span> (polish) {
<a name="l00312"></a>00312                     <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, ADAPTIVE);
<a name="l00313"></a>00313                     <a class="code" href="a00100.html#g0a6fc89686a91e3c8c99923a7b0247e3">evalAndDeflate</a>(p, z, qp);
<a name="l00314"></a>00314                 }
<a name="l00315"></a>00315                 zeros.push_back(z);
<a name="l00316"></a>00316                 *(ComplexPolynomial*)&amp;p = qp;
<a name="l00317"></a>00317             }
<a name="l00318"></a>00318         }
<a name="l00319"></a>00319         <span class="keywordflow">else</span> 
<a name="l00320"></a>00320             <span class="keywordflow">break</span>;
<a name="l00321"></a>00321     }
<a name="l00322"></a>00322     <a class="code" href="a00102.html#g3536daa88b2578a96cd84d2680ac9d66">sortZeros</a>(zeros);
<a name="l00323"></a>00323     <span class="keywordflow">return</span> (<span class="keywordtype">int</span>)zeros.size();
<a name="l00324"></a>00324 }
<a name="l00325"></a>00325 
<a name="l00326"></a>00326 <span class="preprocessor">#endif //JENKINSTRAUB_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>