<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Cohen Acceleration</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Math Toolkit 3.0.0">
<link rel="up" href="../internals.html" title="Internal tools">
<link rel="prev" href="recurrence.html" title="Tools For 3-Term Recurrence Relations">
<link rel="next" href="tuples.html" title="Tuples">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="recurrence.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="tuples.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.internals.cohen_acceleration"></a><a class="link" href="cohen_acceleration.html" title="Cohen Acceleration">Cohen Acceleration</a>
</h3></div></div></div>
<h5>
<a name="math_toolkit.internals.cohen_acceleration.h0"></a>
<span class="phrase"><a name="math_toolkit.internals.cohen_acceleration.synopsis"></a></span><a class="link" href="cohen_acceleration.html#math_toolkit.internals.cohen_acceleration.synopsis">Synopsis</a>
</h5>
<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">cohen_acceleration</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
<span class="keyword">namespace</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="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">G</span><span class="special">></span>
<span class="keyword">auto</span> <span class="identifier">cohen_acceleration</span><span class="special">(</span><span class="identifier">G</span><span class="special">&</span> <span class="identifier">generator</span><span class="special">,</span> <span class="identifier">int64_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="special">-</span><span class="number">1</span><span class="special">);</span>
<span class="special">}</span> <span class="comment">// namespaces</span>
</pre>
<p>
The function <code class="computeroutput"><span class="identifier">cohen_acceleration</span></code>
rapidly computes the limiting value of an alternating series via a technique
developed by <a href="https://www.johndcook.com/blog/2020/08/06/cohen-acceleration/" target="_top">Henri
Cohen et al</a>. To compute
</p>
<p>
<span class="inlinemediaobject"><object type="image/svg+xml" data="../../../equations/alternating_series.svg" width="88" height="51"></object></span>
</p>
<p>
we first define a callable that produces <span class="emphasis"><em>a</em></span><sub><span class="emphasis"><em>k</em></span></sub> on
the kth call. For example, suppose we wish to compute
</p>
<p>
<span class="inlinemediaobject"><object type="image/svg+xml" data="../../../equations/zeta_related_alternating.svg" width="168" height="51"></object></span>
</p>
<p>
First, we need to define a callable which returns the requisite terms:
</p>
<pre class="programlisting"><span class="keyword">template</span><span class="special"><</span><span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">></span>
<span class="keyword">class</span> <span class="identifier">G</span> <span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
<span class="identifier">G</span><span class="special">(){</span>
<span class="identifier">k_</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span>
<span class="special">}</span>
<span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()()</span> <span class="special">{</span>
<span class="identifier">k_</span> <span class="special">+=</span> <span class="number">1</span><span class="special">;</span>
<span class="keyword">return</span> <span class="number">1</span><span class="special">/(</span><span class="identifier">k_</span><span class="special">*</span><span class="identifier">k_</span><span class="special">);</span>
<span class="special">}</span>
<span class="keyword">private</span><span class="special">:</span>
<span class="identifier">Real</span> <span class="identifier">k_</span><span class="special">;</span>
<span class="special">};</span>
</pre>
<p>
Then we pass this into the <code class="computeroutput"><span class="identifier">cohen_acceleration</span></code>
function:
</p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">gen</span> <span class="special">=</span> <span class="identifier">G</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
<span class="keyword">double</span> <span class="identifier">computed</span> <span class="special">=</span> <span class="identifier">cohen_acceleration</span><span class="special">(</span><span class="identifier">gen</span><span class="special">);</span>
</pre>
<p>
See <code class="computeroutput"><span class="identifier">cohen_acceleration</span><span class="special">.</span><span class="identifier">cpp</span></code> in the <code class="computeroutput"><span class="identifier">examples</span></code>
directory for more.
</p>
<p>
The number of terms consumed is computed from the error model
</p>
<p>
<span class="inlinemediaobject"><object type="image/svg+xml" data="../../../equations/cohen_acceleration_error_model.svg" width="394" height="47"></object></span>
</p>
<p>
and must be computed <span class="emphasis"><em>a priori</em></span>. If we read the reference
carefully, we notice that this error model is derived under the assumption
that the terms <span class="emphasis"><em>a</em></span><sub><span class="emphasis"><em>k</em></span></sub> are given as the
moments of a positive measure on [0,1]. If this assumption does not hold,
then the number of terms chosen by the method is incorrect. Hence we permit
the user to provide a second argument to specify the number of terms:
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">computed</span> <span class="special">=</span> <span class="identifier">cohen_acceleration</span><span class="special">(</span><span class="identifier">gen</span><span class="special">,</span> <span class="number">5</span><span class="special">);</span>
</pre>
<p>
<span class="emphasis"><em>Nota bene</em></span>: When experimenting with this option, we found
that adding more terms was no guarantee of additional accuracy, and could
not find an example where a user-provided number of terms outperformed the
default. In addition, it is easy to generate intermediates which overflow
if we let <span class="emphasis"><em>n</em></span> grow too large. Hence we recommend only
playing with this parameter to <span class="emphasis"><em>decrease</em></span> the default
number of terms to increase speed.
</p>
<h5>
<a name="math_toolkit.internals.cohen_acceleration.h1"></a>
<span class="phrase"><a name="math_toolkit.internals.cohen_acceleration.performance"></a></span><a class="link" href="cohen_acceleration.html#math_toolkit.internals.cohen_acceleration.performance">Performance</a>
</h5>
<p>
To see that Cohen acceleration is in fact faster than naive summation for
the same level of relative accuracy, we can run the <code class="computeroutput"><span class="identifier">reporting</span><span class="special">/</span><span class="identifier">performance</span><span class="special">/</span><span class="identifier">cohen_acceleration_performance</span><span class="special">.</span><span class="identifier">cpp</span></code> file.
This benchmark computes the alternating Basel series discussed above:
</p>
<pre class="programlisting"><span class="identifier">Running</span> <span class="special">./</span><span class="identifier">reporting</span><span class="special">/</span><span class="identifier">performance</span><span class="special">/</span><span class="identifier">cohen_acceleration_performance</span><span class="special">.</span><span class="identifier">x</span>
<span class="identifier">Run</span> <span class="identifier">on</span> <span class="special">(</span><span class="number">16</span> <span class="identifier">X</span> <span class="number">2300</span> <span class="identifier">MHz</span> <span class="identifier">CPU</span> <span class="identifier">s</span><span class="special">)</span>
<span class="identifier">CPU</span> <span class="identifier">Caches</span><span class="special">:</span>
<span class="identifier">L1</span> <span class="identifier">Data</span> <span class="number">32</span> <span class="identifier">KiB</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
<span class="identifier">L1</span> <span class="identifier">Instruction</span> <span class="number">32</span> <span class="identifier">KiB</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
<span class="identifier">L2</span> <span class="identifier">Unified</span> <span class="number">256</span> <span class="identifier">KiB</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
<span class="identifier">L3</span> <span class="identifier">Unified</span> <span class="number">16384</span> <span class="identifier">KiB</span> <span class="special">(</span><span class="identifier">x1</span><span class="special">)</span>
<span class="identifier">Load</span> <span class="identifier">Average</span><span class="special">:</span> <span class="number">4.13</span><span class="special">,</span> <span class="number">3.71</span><span class="special">,</span> <span class="number">3.30</span>
<span class="special">-----------------------------------------------------------------</span>
<span class="identifier">Benchmark</span> <span class="identifier">Time</span>
<span class="special">-----------------------------------------------------------------</span>
<span class="identifier">CohenAcceleration</span><span class="special"><</span><span class="keyword">float</span><span class="special">></span> <span class="number">20.7</span> <span class="identifier">ns</span>
<span class="identifier">CohenAcceleration</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="number">64.6</span> <span class="identifier">ns</span>
<span class="identifier">CohenAcceleration</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">></span> <span class="number">115</span> <span class="identifier">ns</span>
<span class="identifier">NaiveSum</span><span class="special"><</span><span class="keyword">float</span><span class="special">></span> <span class="number">4994</span> <span class="identifier">ns</span>
<span class="identifier">NaiveSum</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="number">112803698</span> <span class="identifier">ns</span>
<span class="identifier">NaiveSum</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">></span> <span class="number">5009564877</span> <span class="identifier">ns</span>
</pre>
<p>
In fact not only does the naive sum take orders of magnitude longer to compute,
it is less accurate as well.
</p>
<h5>
<a name="math_toolkit.internals.cohen_acceleration.h2"></a>
<span class="phrase"><a name="math_toolkit.internals.cohen_acceleration.references"></a></span><a class="link" href="cohen_acceleration.html#math_toolkit.internals.cohen_acceleration.references">References</a>
</h5>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "><li class="listitem">
Cohen, Henri, Fernando Rodriguez Villegas, and Don Zagier. <span class="emphasis"><em>Convergence
acceleration of alternating series.</em></span> Experimental mathematics
9.1 (2000): 3-12.
</li></ul></div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="recurrence.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="tuples.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>