📄 a00106.html
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"><html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1"><title>The Polynomials Templates Library: muller.h Source File</title><link href="doxygen.css" rel="stylesheet" type="text/css"><link href="tabs.css" rel="stylesheet" type="text/css"></head><body><!-- Generated by Doxygen 1.4.5 --><div class="tabs"> <ul> <li><a href="index.html"><span>Main Page</span></a></li> <li><a href="modules.html"><span>Modules</span></a></li> <li><a href="annotated.html"><span>Classes</span></a></li> <li id="current"><a href="files.html"><span>Files</span></a></li> </ul></div><div class="tabs"> <ul> <li><a href="files.html"><span>File List</span></a></li> <li><a href="globals.html"><span>File Members</span></a></li> </ul></div><h1>muller.h</h1><a href="a00095.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="preprocessor">#ifndef MULLER_H</span><a name="l00002"></a>00002 <span class="preprocessor"></span><span class="preprocessor">#define MULLER_H</span><a name="l00003"></a>00003 <span class="preprocessor"></span><a name="l00029"></a>00029 <span class="preprocessor">#include "<a class="code" href="a00097.html">polyzero.h</a>"</span><a name="l00030"></a>00030 <a name="l00054"></a>00054 <span class="keyword">template</span> < <span class="keyword">class</span> T, <span class="keyword">class</span> U ><a name="l00055"></a><a class="code" href="a00102.html#g255f8995a5a65456d0ac4b99a1d60325">00055</a> <span class="keywordtype">int</span> <a class="code" href="a00102.html#g255f8995a5a65456d0ac4b99a1d60325">mullerZeros</a>(<span class="keyword">const</span> Polynomial < T > & P,<a name="l00056"></a>00056 std::vector < std::complex < U > > & zeros,<a name="l00057"></a>00057 <span class="keywordtype">bool</span> polish)<a name="l00058"></a>00058 {<a name="l00059"></a>00059 <a name="l00060"></a>00060 <span class="keyword">typedef</span> U Float;<a name="l00061"></a>00061 <span class="keyword">typedef</span> std::complex < Float > Complex;<a name="l00062"></a>00062 <span class="keyword">typedef</span> Polynomial < Complex > ComplexPolynomial;<a name="l00063"></a>00063 <span class="keyword">typedef</span> Polynomial < Float > FloatPolynomial;<a name="l00064"></a>00064 <a name="l00065"></a>00065 <span class="keyword">const</span> <span class="keywordtype">bool</span> ADAPTIVE = <span class="keyword">false</span>;<a name="l00066"></a>00066 <span class="keyword">static</span> FloatSpecs < Float > fpSpecs;<a name="l00067"></a>00067 <span class="keyword">const</span> <span class="keywordtype">int</span> MAXIT = 128;<a name="l00068"></a>00068 <span class="keyword">static</span> <span class="keyword">const</span> Float SHIFT = acos(Float(-1)) * Float(7 / 180.);<a name="l00069"></a>00069 <span class="keyword">const</span> Float GROW = Float(1.05);<a name="l00070"></a>00070 <a name="l00071"></a>00071 Polynomial < T > p = P, f;<a name="l00072"></a>00072 ComplexPolynomial q, q2;<a name="l00073"></a>00073 FloatPolynomial mp;<a name="l00074"></a>00074 <a name="l00075"></a>00075 Complex z0, zx, z, z1, z2, pz, pz1, pz2, bestZ, zE, shift;<a name="l00076"></a>00076 Complex h, h1, h2, d, d1, d2, ac, sq, b, D, E, x, sm, df, px;<a name="l00077"></a>00077 Float e, a = 0, da = 1, r, mpz, mpx, g;<a name="l00078"></a>00078 <span class="keywordtype">int</span> i, j = 0, conv = 0, n;<a name="l00079"></a>00079 <span class="keywordtype">bool</span> started;<a name="l00080"></a>00080 <a name="l00081"></a>00081 zeros.clear();<a name="l00082"></a>00082 <a name="l00083"></a>00083 <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mp);<a name="l00084"></a>00084 <a class="code" href="a00099.html#g663581aef6853ed86df03a35505c9ebd">scalePoly</a>(p, mp);<a name="l00085"></a>00085 f = p;<a name="l00086"></a>00086 i = <a class="code" href="a00102.html#g206394af98d861f1a81bea25ae6d84da">removeNullZeros</a>(p);<a name="l00087"></a>00087 zeros.assign(i, Complex(0));<a name="l00088"></a>00088 <a name="l00089"></a>00089 <span class="keywordflow">while</span> (p.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>() > 0) {<a name="l00090"></a>00090 <a name="l00091"></a>00091 <span class="keywordflow">if</span> (p.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>() == 1) {<a name="l00092"></a>00092 z = <a class="code" href="a00102.html#g6e9fe1cec8c7ec334d9e21f31b602bde">solveDegree1</a>(p[1], p[0]);<a name="l00093"></a>00093 zeros.push_back(z);<a name="l00094"></a>00094 <span class="keywordflow">break</span>;<a name="l00095"></a>00095 }<a name="l00096"></a>00096 <span class="keywordflow">if</span> (p.degree() == 2) {<a name="l00097"></a>00097 <a class="code" href="a00102.html#gb88e92200690d2e3f8b0c785c8fa42d3">solveDegree2</a>(p[2], p[1], p[0], z, z1);<a name="l00098"></a>00098 <span class="keywordflow">if</span> (polish) {<a name="l00099"></a>00099 <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, ADAPTIVE);<a name="l00100"></a>00100 <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z1, pz, mpz, ADAPTIVE);<a name="l00101"></a>00101 }<a name="l00102"></a>00102 zeros.push_back(z);<a name="l00103"></a>00103 zeros.push_back(z1);<a name="l00104"></a>00104 <span class="keywordflow">break</span>;<a name="l00105"></a>00105 }<a name="l00106"></a>00106 <a name="l00107"></a>00107 <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mp);<a name="l00108"></a>00108 mp[0] = -mp[0];<a name="l00109"></a>00109 Float c = <a class="code" href="a00099.html#g4c00b88d043d1f70cdcb4d2cec66b55e">cauchyLowerBound</a>(mp, Float(0));<a name="l00110"></a>00110 r = c;<a name="l00111"></a>00111 g = <a class="code" href="a00099.html#ga911d353ac1900e9da22160a11d8b350">zerosGeometricMean</a>(p);<a name="l00112"></a>00112 g = (g - c) / MAXIT;<a name="l00113"></a>00113 <span class="keywordflow">for</span> (i = 1; i <= MAXIT; ++i) {<a name="l00114"></a>00114 <a name="l00115"></a>00115 shift = exp(Complex(0, a));<a name="l00116"></a>00116 z2 = Complex(0, -r) * shift;<a name="l00117"></a>00117 z1 = 0;<a name="l00118"></a>00118 z = (r) * shift;<a name="l00119"></a>00119 pz = <a class="code" href="a00100.html#gc5161c8ecd75b0108d0ea60c92d30d21">eval</a>(p, z);<a name="l00120"></a>00120 pz1 = <a class="code" href="a00100.html#gc5161c8ecd75b0108d0ea60c92d30d21">eval</a>(p, z1);<a name="l00121"></a>00121 pz2 = <a class="code" href="a00100.html#gc5161c8ecd75b0108d0ea60c92d30d21">eval</a>(p, z2);<a name="l00122"></a>00122 <span class="keywordflow">if</span> (abs(pz2) < abs(pz)) {<a name="l00123"></a>00123 std::swap(pz, pz2);<a name="l00124"></a>00124 std::swap(z, z2);<a name="l00125"></a>00125 }<a name="l00126"></a>00126 <a name="l00127"></a>00127 a += SHIFT;<a name="l00128"></a>00128 r += g;<a name="l00129"></a>00129 <a name="l00130"></a>00130 e = abs(pz) * (1 + fpSpecs.MRE), mpx, mpz;<a name="l00131"></a>00131 started = <span class="keyword">false</span>;<a name="l00132"></a>00132 <a name="l00133"></a>00133 h1 = z1 - z2;<a name="l00134"></a>00134 h2 = z - z1;<a name="l00135"></a>00135 <a name="l00136"></a>00136 <span class="keywordflow">if</span> (h1 == Complex(0) && h2 == Complex(0) || (h1 + h2) == Complex(0))<a name="l00137"></a>00137 <span class="keywordflow">continue</span>;<a name="l00138"></a>00138 d1 = (pz1 - pz2) / h1;<a name="l00139"></a>00139 d2 = (pz - pz1) / h2;<a name="l00140"></a>00140 d = (d2 - d1) / (h2 + h1);<a name="l00141"></a>00141 <a name="l00142"></a>00142 px = pz;<a name="l00143"></a>00143 mpx = mpz = abs(pz);<a name="l00144"></a>00144 <a name="l00145"></a>00145 j = 0;<a name="l00146"></a>00146 conv = 0;<a name="l00147"></a>00147 <span class="keywordflow">while</span> (1) {<a name="l00148"></a>00148 <span class="keywordflow">if</span> (++j > 3)<a name="l00149"></a>00149 started = <span class="keyword">true</span>;<a name="l00150"></a>00150 <a name="l00151"></a>00151 <span class="keywordflow">if</span> (!started || conv == 0) {<a name="l00152"></a>00152 b = d2 + (h2 * d);<a name="l00153"></a>00153 ac = pz * d;<a name="l00154"></a>00154 ac = ac + ac;<a name="l00155"></a>00155 ac = ac + ac;<a name="l00156"></a>00156 D = sqrt((b * b) - ac);<a name="l00157"></a>00157 <a name="l00158"></a>00158 sm = b + D;<a name="l00159"></a>00159 df = b - D;<a name="l00160"></a>00160 E = (abs(df) > abs(sm)) ? df : sm;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -