📄 a00105.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: laguerre.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>laguerre.h</h1><a href="a00094.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="preprocessor">#ifndef LAGUERRE_H</span><a name="l00002"></a>00002 <span class="preprocessor"></span><span class="preprocessor">#define LAGUERRE_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="l00061"></a>00061 <span class="keyword">template</span><<span class="keyword">class</span> T, <span class="keyword">class</span> U><a name="l00062"></a><a class="code" href="a00102.html#g8f609e24d651d3dcf78a90ebf4a92380">00062</a> <span class="keywordtype">int</span> <a class="code" href="a00102.html#g8f609e24d651d3dcf78a90ebf4a92380">laguerreZeros</a>(<span class="keyword">const</span> Polynomial< T >& P, std::vector< std::complex < U > >& zeros, <span class="keywordtype">bool</span> polish) {<a name="l00063"></a>00063 <a name="l00064"></a>00064 <span class="keyword">typedef</span> U Float;<a name="l00065"></a>00065 <span class="keyword">typedef</span> std::complex<Float> Complex;<a name="l00066"></a>00066 <span class="keyword">typedef</span> Polynomial<Complex> ComplexPolynomial;<a name="l00067"></a>00067 <span class="keyword">typedef</span> Polynomial<Float> FloatPolynomial;<a name="l00068"></a>00068 <a name="l00069"></a>00069 <span class="comment">// local scope functions </span><a name="l00070"></a>00070 <span class="keyword">class </span>local {<a name="l00071"></a>00071 <span class="keyword">public</span>:<a name="l00072"></a>00072 <span class="keyword">static</span> Complex fejer(Float m, <span class="keyword">const</span> Complex& pz, <span class="keyword">const</span> Complex& ppz, <span class="keyword">const</span> Complex& pppz) {<a name="l00073"></a>00073 Complex t1, t2, t3, z1, z2;<a name="l00074"></a>00074 t1 = pppz / (m * (m - 1));<a name="l00075"></a>00075 t2 = Float(2) * ppz / m;<a name="l00076"></a>00076 t3 = pz;<a name="l00077"></a>00077 <a class="code" href="a00102.html#gb88e92200690d2e3f8b0c785c8fa42d3">solveDegree2</a>(t1, t2, t3, z1, z2);<a name="l00078"></a>00078 <span class="keywordflow">return</span> z1;<a name="l00079"></a>00079 }<a name="l00080"></a>00080 <a name="l00081"></a>00081 <span class="keyword">static</span> Complex laguerreStep(Float m, <span class="keyword">const</span> Complex& zs, <span class="keyword">const</span> Complex& pz, <span class="keyword">const</span> Complex& ppz) {<a name="l00082"></a>00082 Complex t;<a name="l00083"></a>00083 t = ((m - 2) * ppz) / (m * pz);<a name="l00084"></a>00084 t = (zs * t) + (m - 1);<a name="l00085"></a>00085 <span class="keywordflow">return</span> zs / t;<a name="l00086"></a>00086 }<a name="l00087"></a>00087 };<a name="l00088"></a>00088 <a name="l00089"></a>00089 <span class="keyword">static</span> FloatSpecs<Float> fpSpecs;<a name="l00090"></a>00090 <span class="keyword">const</span> Float SMALL = Float(1.e-3);<a name="l00091"></a>00091 <span class="keyword">const</span> Float BIG_ONE = Float(1.0001);<a name="l00092"></a>00092 <span class="keyword">const</span> Float SML_ONE = Float(0.9999);<a name="l00093"></a>00093 <span class="keyword">const</span> Float RCONST = Float(1.445);<a name="l00094"></a>00094 <span class="keyword">const</span> Float ONEPQT = Float(1.25);<a name="l00095"></a>00095 <span class="keyword">const</span> Float GAMA = Float(1);<a name="l00096"></a>00096 <span class="keyword">const</span> Float THETA = Float(2);<a name="l00097"></a>00097 <a name="l00098"></a>00098 Polynomial<T> p = P;<a name="l00099"></a>00099 ComplexPolynomial q, q2;<a name="l00100"></a>00100 FloatPolynomial mP, d, p2, rem;<a name="l00101"></a>00101 <span class="keywordtype">int</span> startd, iter, spiral, n;<a name="l00102"></a>00102 Complex spir, temp, r;<a name="l00103"></a>00103 Float sqm, m, ratio;<a name="l00104"></a>00104 <a name="l00105"></a>00105 Complex zs, l0, step, l, z0, z, pz, ppz, pppz, pz0;<a name="l00106"></a>00106 Float c, f, g, w, ml0, ub, lb, mstep, ml, mpz, mpz0;<a name="l00107"></a>00107 <a name="l00108"></a>00108 d.resize(3);<a name="l00109"></a>00109 <a name="l00110"></a>00110 <a class="code" href="a00099.html#g1ca5c1a34fcdcadb083a25ecf0e1a995">modulus</a>(p, mP);<a name="l00111"></a>00111 <a class="code" href="a00099.html#g663581aef6853ed86df03a35505c9ebd">scalePoly</a>(p, mP);<a name="l00112"></a>00112 mP[0] = -mP[0];<a name="l00113"></a>00113 <a name="l00114"></a>00114 n = <a class="code" href="a00102.html#g206394af98d861f1a81bea25ae6d84da">removeNullZeros</a>(p);<a name="l00115"></a>00115 zeros.assign(n, Complex(0));<a name="l00116"></a>00116 <a name="l00117"></a>00117 <span class="keywordflow">while</span> (p.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>() > 0) {<a name="l00118"></a>00118 <a name="l00119"></a>00119 m = Float(p.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>());<a name="l00120"></a>00120 <a name="l00121"></a>00121 <span class="keywordflow">if</span> (p.<a class="code" href="a00090.html#b433e0edb11b32a205a9d07b89aefa8f">degree</a>() == 1) {<a name="l00122"></a>00122 z = <a class="code" href="a00102.html#g6e9fe1cec8c7ec334d9e21f31b602bde">solveDegree1</a>(p[1], p[0]);<a name="l00123"></a>00123 zeros.push_back(z);<a name="l00124"></a>00124 <span class="keywordflow">break</span>;<a name="l00125"></a>00125 }<a name="l00126"></a>00126 <span class="keywordflow">if</span> (p.degree() == 2) {<a name="l00127"></a>00127 <a class="code" href="a00102.html#gb88e92200690d2e3f8b0c785c8fa42d3">solveDegree2</a>(p[2], p[1], p[0], z0, z);<a name="l00128"></a>00128 <span class="keywordflow">if</span> (polish) {<a name="l00129"></a>00129 <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z, pz, mpz, <span class="keyword">false</span>);<a name="l00130"></a>00130 <a class="code" href="a00102.html#g7528de40f73b280a95202bc5b8947424">newtonZero</a>(P, z0, pz, mpz, <span class="keyword">false</span>);<a name="l00131"></a>00131 }<a name="l00132"></a>00132 zeros.push_back(z0);<a name="l00133"></a>00133 zeros.push_back(z);<a name="l00134"></a>00134 <span class="keywordflow">break</span>;<a name="l00135"></a>00135 }<a name="l00136"></a>00136 <a name="l00137"></a>00137 <span class="comment">// get an initial estimate z, p(z)</span><a name="l00138"></a>00138 z = Complex(0);<a name="l00139"></a>00139 spir = Complex(-ONEPQT / 10, 1);<a name="l00140"></a>00140 <span class="keywordflow">do</span> {<a name="l00141"></a>00141 pz = <a class="code" href="a00100.html#gc143079d7e07a2b902831ec539e45fc6">evalAndDerive</a>(p, z, ppz, pppz);<a name="l00142"></a>00142 <span class="comment">// compute zs, fejer bound</span><a name="l00143"></a>00143 zs = local::fejer(m, pz, ppz, pppz);<a name="l00144"></a>00144 f = abs(zs);<a name="l00145"></a>00145 <a name="l00146"></a>00146 <span class="comment">// diro = l0</span><a name="l00147"></a>00147 l0 = local::laguerreStep(m, zs, p[0], p[1]);<a name="l00148"></a>00148 ml0 = abs(l0);<a name="l00149"></a>00149 <a name="l00150"></a>00150 <span class="keywordflow">if</span> (_isnan(ml0)) {<a name="l00151"></a>00151 <span class="keywordflow">if</span> (z == Complex(0))<a name="l00152"></a>00152 z = -p[0];<a name="l00153"></a>00153 <span class="keywordflow">else</span><a name="l00154"></a>00154 z = spir * z;<a name="l00155"></a>00155 }<a name="l00156"></a>00156 <span class="keywordflow">else</span><a name="l00157"></a>00157 <span class="keywordflow">break</span>;<a name="l00158"></a>00158 } <span class="keywordflow">while</span> (1);<a name="l00159"></a>00159 <a name="l00160"></a>00160 sqm = sqrt(m);<a name="l00161"></a>00161 g = <a class="code" href="a00099.html#ga911d353ac1900e9da22160a11d8b350">zerosGeometricMean</a>(p);<a name="l00162"></a>00162 <a name="l00163"></a>00163 <span class="comment">// Laguerre bound w = sqrt(m).|step|;</span><a name="l00164"></a>00164 w = sqm * ml0;<a name="l00165"></a>00165 <a name="l00166"></a>00166 <span class="comment">// get cauchy bound</span><a name="l00167"></a>00167 ub = std::min(g, BIG_ONE * std::min(f, w));<a name="l00168"></a>00168 c = <a class="code" href="a00099.html#g4c00b88d043d1f70cdcb4d2cec66b55e">cauchyLowerBound</a>(mP, ub);<a name="l00169"></a>00169 <a name="l00170"></a>00170 <span class="comment">// compute upper and lower bound of smallest zero</span><a name="l00171"></a>00171 ub = std::min(RCONST * m * c, ub);<a name="l00172"></a>00172 lb = SML_ONE * c;<a name="l00173"></a>00173 <a name="l00174"></a>00174 <span class="comment">// init first pass</span><a name="l00175"></a>00175 f = ub;<a name="l00176"></a>00176 g = ub;<a name="l00177"></a>00177 step = l0;<a name="l00178"></a>00178 mstep = ml0;<a name="l00179"></a>00179 ratio = mstep / g;<a name="l00180"></a>00180 mpz = abs(pz);<a name="l00181"></a>00181 mpz0 = mpz;<a name="l00182"></a>00182 spiral = 0;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -