1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Series Evaluation</title> 5<link rel="stylesheet" href="../../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../internals.html" title="Internal tools"> 9<link rel="prev" href="../internals.html" title="Internal tools"> 10<link rel="next" href="cf.html" title="Continued Fraction Evaluation"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="../internals.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../internals.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="cf.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.internals.series_evaluation"></a><a class="link" href="series_evaluation.html" title="Series Evaluation">Series Evaluation</a> 28</h3></div></div></div> 29<h5> 30<a name="math_toolkit.internals.series_evaluation.h0"></a> 31 <span class="phrase"><a name="math_toolkit.internals.series_evaluation.synopsis"></a></span><a class="link" href="series_evaluation.html#math_toolkit.internals.series_evaluation.synopsis">Synopsis</a> 32 </h5> 33<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">series</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34</pre> 35<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">tools</span><span class="special">{</span> 36 37<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">U</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">V</span><span class="special">></span> 38<span class="keyword">inline</span> <span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">U</span><span class="special">&</span> <span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_terms</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">V</span><span class="special">&</span> <span class="identifier">init_value</span><span class="special">);</span> 39 40<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">U</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">V</span><span class="special">></span> 41<span class="keyword">inline</span> <span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">const</span> <span class="identifier">U</span><span class="special">&</span> <span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_terms</span><span class="special">);</span> 42 43<span class="comment">//</span> 44<span class="comment">// The following interfaces are now deprecated:</span> 45<span class="comment">// </span> 46<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">></span> 47<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">);</span> 48 49<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">></span> 50<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_terms</span><span class="special">);</span> 51 52<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">U</span><span class="special">></span> 53<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">U</span> <span class="identifier">init_value</span><span class="special">);</span> 54 55<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">U</span><span class="special">></span> 56<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_terms</span><span class="special">,</span> <span class="identifier">U</span> <span class="identifier">init_value</span><span class="special">);</span> 57 58<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">></span> 59<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">kahan_sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">);</span> 60 61<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Functor</span><span class="special">></span> 62<span class="keyword">typename</span> <span class="identifier">Functor</span><span class="special">::</span><span class="identifier">result_type</span> <span class="identifier">kahan_sum_series</span><span class="special">(</span><span class="identifier">Functor</span><span class="special">&</span> <span class="identifier">func</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_terms</span><span class="special">);</span> 63 64<span class="special">}}}</span> <span class="comment">// namespaces</span> 65</pre> 66<h5> 67<a name="math_toolkit.internals.series_evaluation.h1"></a> 68 <span class="phrase"><a name="math_toolkit.internals.series_evaluation.description"></a></span><a class="link" href="series_evaluation.html#math_toolkit.internals.series_evaluation.description">Description</a> 69 </h5> 70<p> 71 These algorithms are intended for the <a href="http://en.wikipedia.org/wiki/Series_%28mathematics%29" target="_top">summation 72 of infinite series</a>. 73 </p> 74<p> 75 Each of the algorithms takes a nullary-function object as the first argument: 76 the function object will be repeatedly invoked to pull successive terms from 77 the series being summed. 78 </p> 79<p> 80 The second argument is the precision required, summation will stop when the 81 next term is less than <span class="emphasis"><em>tolerance</em></span> times the result. The 82 deprecated versions of <code class="computeroutput"><span class="identifier">sum_series</span></code> 83 take an integer number of bits here - internally they just convert this to 84 a tolerance and forward the call. 85 </p> 86<p> 87 The third argument <span class="emphasis"><em>max_terms</em></span> sets an upper limit on 88 the number of terms of the series to evaluate. In addition, on exit the function 89 will set <span class="emphasis"><em>max_terms</em></span> to the actual number of terms of 90 the series that were evaluated: this is particularly useful for profiling 91 the convergence properties of a new series. 92 </p> 93<p> 94 The final optional argument <span class="emphasis"><em>init_value</em></span> is the initial 95 value of the sum to which the terms of the series should be added. This is 96 useful in two situations: 97 </p> 98<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 99<li class="listitem"> 100 Where the first value of the series has a different formula to successive 101 terms. In this case the first value in the series can be passed as the 102 last argument and the logic of the function object can then be simplified 103 to return subsequent terms. 104 </li> 105<li class="listitem"> 106 Where the series is being added (or subtracted) from some other value: 107 termination of the series will likely occur much more rapidly if that 108 other value is passed as the last argument. For example, there are several 109 functions that can be expressed as <span class="emphasis"><em>1 - S(z)</em></span> where 110 S(z) is an infinite series. In this case, pass -1 as the last argument 111 and then negate the result of the summation to get the result of <span class="emphasis"><em>1 112 - S(z)</em></span>. 113 </li> 114</ul></div> 115<p> 116 The two <span class="emphasis"><em>kahan_sum_series</em></span> variants of these algorithms 117 maintain a carry term that corrects for roundoff error during summation. 118 They are inspired by the <a href="http://en.wikipedia.org/wiki/Kahan_Summation_Algorithm" target="_top"><span class="emphasis"><em>Kahan 119 Summation Formula</em></span></a> that appears in <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html" target="_top">What 120 Every Computer Scientist Should Know About Floating-Point Arithmetic</a>. 121 However, it should be pointed out that there are very few series that require 122 summation in this way. 123 </p> 124<h5> 125<a name="math_toolkit.internals.series_evaluation.h2"></a> 126 <span class="phrase"><a name="math_toolkit.internals.series_evaluation.examples"></a></span><a class="link" href="series_evaluation.html#math_toolkit.internals.series_evaluation.examples">Examples</a> 127 </h5> 128<p> 129 These examples are all in <a href="../../../../example/series.cpp" target="_top">../../example/series.cpp</a> 130 </p> 131<p> 132 Let's suppose we want to implement <span class="emphasis"><em>log(1+x)</em></span> via its 133 infinite series, 134 </p> 135<div class="blockquote"><blockquote class="blockquote"><p> 136 <span class="inlinemediaobject"><img src="../../../equations/log1pseries.svg"></span> 137 138 </p></blockquote></div> 139<p> 140 We begin by writing a small function object to return successive terms of 141 the series: 142 </p> 143<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 144<span class="keyword">struct</span> <span class="identifier">log1p_series</span> 145<span class="special">{</span> 146 <span class="comment">// we must define a result_type typedef:</span> 147 <span class="keyword">typedef</span> <span class="identifier">T</span> <span class="identifier">result_type</span><span class="special">;</span> 148 149 <span class="identifier">log1p_series</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 150 <span class="special">:</span> <span class="identifier">k</span><span class="special">(</span><span class="number">0</span><span class="special">),</span> <span class="identifier">m_mult</span><span class="special">(-</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">m_prod</span><span class="special">(-</span><span class="number">1</span><span class="special">)</span> <span class="special">{}</span> 151 152 <span class="identifier">T</span> <span class="keyword">operator</span><span class="special">()()</span> 153 <span class="special">{</span> 154 <span class="comment">// This is the function operator invoked by the summation</span> 155 <span class="comment">// algorithm, the first call to this operator should return</span> 156 <span class="comment">// the first term of the series, the second call the second</span> 157 <span class="comment">// term and so on.</span> 158 <span class="identifier">m_prod</span> <span class="special">*=</span> <span class="identifier">m_mult</span><span class="special">;</span> 159 <span class="keyword">return</span> <span class="identifier">m_prod</span> <span class="special">/</span> <span class="special">++</span><span class="identifier">k</span><span class="special">;</span> 160 <span class="special">}</span> 161 162<span class="keyword">private</span><span class="special">:</span> 163 <span class="keyword">int</span> <span class="identifier">k</span><span class="special">;</span> 164 <span class="keyword">const</span> <span class="identifier">T</span> <span class="identifier">m_mult</span><span class="special">;</span> 165 <span class="identifier">T</span> <span class="identifier">m_prod</span><span class="special">;</span> 166<span class="special">};</span> 167</pre> 168<p> 169 Implementing log(1+x) is now fairly trivial: 170 </p> 171<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 172<span class="identifier">T</span> <span class="identifier">log1p</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 173<span class="special">{</span> 174 <span class="comment">// We really should add some error checking on x here!</span> 175 <span class="identifier">BOOST_ASSERT</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">fabs</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><</span> <span class="number">1</span><span class="special">);</span> 176 177 <span class="comment">// Construct the series functor:</span> 178 <span class="identifier">log1p_series</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">s</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> 179 <span class="comment">// Set a limit on how many iterations we permit:</span> 180 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">max_iter</span> <span class="special">=</span> <span class="number">1000</span><span class="special">;</span> 181 <span class="comment">// Add it up, with enough precision for full machine precision:</span> 182 <span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">(),</span> <span class="identifier">max_iter</span><span class="special">);</span> 183<span class="special">}</span> 184</pre> 185<p> 186 We can almost use the code above for complex numbers as well - unfortunately 187 we need a slightly different definition for epsilon, and within the functor, 188 mixed complex and integer arithmetic is sadly not supported (as of C++17), 189 so we need to cast out integers to floats: 190 </p> 191<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 192<span class="keyword">struct</span> <span class="identifier">log1p_series</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="special">></span> 193<span class="special">{</span> 194 <span class="comment">// we must define a result_type typedef:</span> 195 <span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">result_type</span><span class="special">;</span> 196 197 <span class="identifier">log1p_series</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">x</span><span class="special">)</span> 198 <span class="special">:</span> <span class="identifier">k</span><span class="special">(</span><span class="number">0</span><span class="special">),</span> <span class="identifier">m_mult</span><span class="special">(-</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">m_prod</span><span class="special">(-</span><span class="number">1</span><span class="special">)</span> <span class="special">{}</span> 199 200 <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="keyword">operator</span><span class="special">()()</span> 201 <span class="special">{</span> 202 <span class="comment">// This is the function operator invoked by the summation</span> 203 <span class="comment">// algorithm, the first call to this operator should return</span> 204 <span class="comment">// the first term of the series, the second call the second</span> 205 <span class="comment">// term and so on.</span> 206 <span class="identifier">m_prod</span> <span class="special">*=</span> <span class="identifier">m_mult</span><span class="special">;</span> 207 <span class="keyword">return</span> <span class="identifier">m_prod</span> <span class="special">/</span> <span class="identifier">T</span><span class="special">(++</span><span class="identifier">k</span><span class="special">);</span> 208 <span class="special">}</span> 209 210<span class="keyword">private</span><span class="special">:</span> 211 <span class="keyword">int</span> <span class="identifier">k</span><span class="special">;</span> 212 <span class="keyword">const</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">m_mult</span><span class="special">;</span> 213 <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">m_prod</span><span class="special">;</span> 214<span class="special">};</span> 215 216 217<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 218<span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">log1p</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">x</span><span class="special">)</span> 219<span class="special">{</span> 220 <span class="comment">// We really should add some error checking on x here!</span> 221 <span class="identifier">BOOST_ASSERT</span><span class="special">(</span><span class="identifier">abs</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><</span> <span class="number">1</span><span class="special">);</span> 222 223 <span class="comment">// Construct the series functor:</span> 224 <span class="identifier">log1p_series</span><span class="special"><</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="special">></span> <span class="identifier">s</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> 225 <span class="comment">// Set a limit on how many iterations we permit:</span> 226 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">max_iter</span> <span class="special">=</span> <span class="number">1000</span><span class="special">;</span> 227 <span class="comment">// Add it up, with enough precision for full machine precision:</span> 228 <span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">sum_series</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">complex</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()),</span> <span class="identifier">max_iter</span><span class="special">);</span> 229<span class="special">}</span> 230</pre> 231<p> 232 Of course with a few traits classes and a bit of meta-programming we could 233 fold these two implementations into one, but that's beyond the scope of these 234 examples. 235 </p> 236</div> 237<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 238<td align="left"></td> 239<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 240 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 241 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 242 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 243 Daryle Walker and Xiaogang Zhang<p> 244 Distributed under the Boost Software License, Version 1.0. (See accompanying 245 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 246 </p> 247</div></td> 248</tr></table> 249<hr> 250<div class="spirit-nav"> 251<a accesskey="p" href="../internals.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../internals.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="cf.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 252</div> 253</body> 254</html> 255
