1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Cardinal Cubic B-spline interpolation</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="../interpolation.html" title="Chapter 12. Interpolation"> 9<link rel="prev" href="../interpolation.html" title="Chapter 12. Interpolation"> 10<link rel="next" href="cardinal_quadratic_b.html" title="Cardinal Quadratic B-spline interpolation"> 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="../interpolation.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../interpolation.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="cardinal_quadratic_b.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.cardinal_cubic_b"></a><a class="link" href="cardinal_cubic_b.html" title="Cardinal Cubic B-spline interpolation">Cardinal Cubic B-spline 28 interpolation</a> 29</h2></div></div></div> 30<h4> 31<a name="math_toolkit.cardinal_cubic_b.h0"></a> 32 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.synopsis"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.synopsis">Synopsis</a> 33 </h4> 34<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">interpolators</span><span class="special">/</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 35</pre> 36<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">interpolators</span> <span class="special">{</span> 37 38 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span> 39 <span class="keyword">class</span> <span class="identifier">cardinal_cubic_b_spline</span> 40 <span class="special">{</span> 41 <span class="keyword">public</span><span class="special">:</span> 42 43 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">BidiIterator</span><span class="special">></span> 44 <span class="identifier">cardinal_cubic_b_spline</span><span class="special">(</span><span class="identifier">BidiIterator</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">BidiIterator</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">left_endpoint</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">step_size</span><span class="special">,</span> 45 <span class="identifier">Real</span> <span class="identifier">left_endpoint_derivative</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">Real</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">(),</span> 46 <span class="identifier">Real</span> <span class="identifier">right_endpoint_derivative</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">Real</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">());</span> 47 <span class="identifier">cardinal_cubic_b_spline</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">Real</span><span class="special">*</span> <span class="keyword">const</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">size_t</span> <span class="identifier">length</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">left_endpoint</span><span class="special">,</span> <span class="identifier">Real</span> <span class="identifier">step_size</span><span class="special">,</span> 48 <span class="identifier">Real</span> <span class="identifier">left_endpoint_derivative</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">Real</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">(),</span> 49 <span class="identifier">Real</span> <span class="identifier">right_endpoint_derivative</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">Real</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">());</span> 50 51 <span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 52 53 <span class="identifier">Real</span> <span class="identifier">prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 54 55 <span class="identifier">Real</span> <span class="identifier">double_prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 56 <span class="special">};</span> 57 58<span class="special">}}}</span> <span class="comment">// namespaces</span> 59</pre> 60<h4> 61<a name="math_toolkit.cardinal_cubic_b.h1"></a> 62 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.cardinal_cubic_b_spline_interpol"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.cardinal_cubic_b_spline_interpol">Cardinal 63 Cubic B-Spline Interpolation</a> 64 </h4> 65<p> 66 The cardinal cubic <span class="emphasis"><em>B</em></span>-spline class provided by Boost allows 67 fast and accurate interpolation of a function which is known at equally spaced 68 points. The cubic <span class="emphasis"><em>B</em></span>-spline interpolation is numerically 69 stable as it uses compactly supported basis functions constructed via iterative 70 convolution. This is to be contrasted to one-sided power function cubic spline 71 interpolation which is ill-conditioned as the global support of cubic polynomials 72 causes small changes far from the evaluation point to exert a large influence 73 on the calculated value. 74 </p> 75<p> 76 There are many use cases for interpolating a function at equally spaced points. 77 One particularly important example is solving ODEs whose coefficients depend 78 on data determined from experiment or numerical simulation. Since most ODE 79 steppers are adaptive, they must be able to sample the coefficients at arbitrary 80 points; not just at the points we know the values of our function. 81 </p> 82<p> 83 The first two arguments to the constructor are either: 84 </p> 85<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 86<li class="listitem"> 87 A pair of bidirectional iterators into the data, or 88 </li> 89<li class="listitem"> 90 A pointer to the data, and a length of the data array. 91 </li> 92</ul></div> 93<p> 94 These are then followed by: 95 </p> 96<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 97<li class="listitem"> 98 The start of the functions domain, 99 </li> 100<li class="listitem"> 101 The step size. 102 </li> 103</ul></div> 104<p> 105 Optionally, you may provide two additional arguments to the constructor, namely 106 the derivative of the function at the left endpoint, and the derivative at 107 the right endpoint. If you do not provide these arguments, they will be estimated 108 using one-sided finite-difference formulas. An example of a valid call to the 109 constructor is 110 </p> 111<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">f</span><span class="special">{</span><span class="number">0.01</span><span class="special">,</span> <span class="special">-</span><span class="number">0.02</span><span class="special">,</span> <span class="number">0.3</span><span class="special">,</span> <span class="number">0.8</span><span class="special">,</span> <span class="number">1.9</span><span class="special">,</span> <span class="special">-</span><span class="number">8.78</span><span class="special">,</span> <span class="special">-</span><span class="number">22.6</span><span class="special">};</span> 112<span class="keyword">double</span> <span class="identifier">t0</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> 113<span class="keyword">double</span> <span class="identifier">h</span> <span class="special">=</span> <span class="number">0.01</span><span class="special">;</span> 114<span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">interpolators</span><span class="special">::</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">spline</span><span class="special">(</span><span class="identifier">f</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">f</span><span class="special">.</span><span class="identifier">end</span><span class="special">(),</span> <span class="identifier">t0</span><span class="special">,</span> <span class="identifier">h</span><span class="special">);</span> 115</pre> 116<p> 117 The endpoints are estimated using a one-sided finite-difference formula. If 118 you know the derivative at the endpoint, you may pass it to the constructor 119 via 120 </p> 121<pre class="programlisting"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">interpolators</span><span class="special">::</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">spline</span><span class="special">(</span><span class="identifier">f</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">f</span><span class="special">.</span><span class="identifier">end</span><span class="special">(),</span> <span class="identifier">t0</span><span class="special">,</span> <span class="identifier">h</span><span class="special">,</span> <span class="identifier">a_prime</span><span class="special">,</span> <span class="identifier">b_prime</span><span class="special">);</span> 122</pre> 123<p> 124 To evaluate the interpolant at a point, we simply use 125 </p> 126<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">spline</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> 127</pre> 128<p> 129 and to evaluate the derivative of the interpolant we use 130 </p> 131<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">yp</span> <span class="special">=</span> <span class="identifier">spline</span><span class="special">.</span><span class="identifier">prime</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> 132</pre> 133<p> 134 Be aware that the accuracy guarantees on the derivative of the spline are an 135 order lower than the guarantees on the original function, see <a href="http://www.springer.com/us/book/9780387984087" target="_top">Numerical 136 Analysis, Graduate Texts in Mathematics, 181, Rainer Kress</a> for details. 137 </p> 138<p> 139 The last interesting member is the second derivative, evaluated via 140 </p> 141<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">ypp</span> <span class="special">=</span> <span class="identifier">spline</span><span class="special">.</span><span class="identifier">double_prime</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> 142</pre> 143<p> 144 The basis functions of the spline are cubic polynomials, so the second derivative 145 is simply linear interpolation. But the interpolation is not constrained by 146 data (you don't pass in the second derivatives), and hence is less accurate 147 than would be naively expected from a linear interpolator. The problem is especially 148 pronounced at the boundaries, where the second derivative is totally unconstrained. 149 Use the second derivative of the cubic B-spline interpolator only in desperation. 150 The quintic <span class="emphasis"><em>B</em></span>-spline interpolator is recommended for cases 151 where second derivatives are needed. 152 </p> 153<h4> 154<a name="math_toolkit.cardinal_cubic_b.h2"></a> 155 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.complexity_and_performance"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.complexity_and_performance">Complexity 156 and Performance</a> 157 </h4> 158<p> 159 The call to the constructor requires ��(<span class="emphasis"><em>n</em></span>) operations, where 160 <span class="emphasis"><em>n</em></span> is the number of points to interpolate. Each call the 161 the interpolant is ��(1) (constant time). On the author's Intel Xeon E3-1230, 162 this takes 21ns as long as the vector is small enough to fit in cache. 163 </p> 164<h4> 165<a name="math_toolkit.cardinal_cubic_b.h3"></a> 166 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.accuracy"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.accuracy">Accuracy</a> 167 </h4> 168<p> 169 Let <span class="emphasis"><em>h</em></span> be the stepsize. If <span class="emphasis"><em>f</em></span> is four-times 170 continuously differentiable, then the interpolant is <span class="emphasis"><em>��(h<sup>4</sup>)</em></span> 171 accurate and the derivative is <span class="emphasis"><em>��(h<sup>3</sup>)</em></span> accurate. 172 </p> 173<h4> 174<a name="math_toolkit.cardinal_cubic_b.h4"></a> 175 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.testing"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.testing">Testing</a> 176 </h4> 177<p> 178 Since the interpolant obeys <span class="serif_italic">s(x<sub>j</sub>) = f(x<sub>j</sub>)</span> 179 at all interpolation points, the tests generate random data and evaluate the 180 interpolant at the interpolation points, validating that equality with the 181 data holds. 182 </p> 183<p> 184 In addition, constant, linear, and quadratic functions are interpolated to 185 ensure that the interpolant behaves as expected. 186 </p> 187<h4> 188<a name="math_toolkit.cardinal_cubic_b.h5"></a> 189 <span class="phrase"><a name="math_toolkit.cardinal_cubic_b.example"></a></span><a class="link" href="cardinal_cubic_b.html#math_toolkit.cardinal_cubic_b.example">Example</a> 190 </h4> 191<p> 192 This example demonstrates how to use the cubic b spline interpolator for regularly 193 spaced data. 194 </p> 195<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">interpolators</span><span class="special">/</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 196 197<span class="keyword">int</span> <span class="identifier">main</span><span class="special">()</span> 198<span class="special">{</span> 199 <span class="comment">// We begin with an array of samples:</span> 200 <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 201 <span class="comment">// And decide on a stepsize:</span> 202 <span class="keyword">double</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">0.01</span><span class="special">;</span> 203 204 <span class="comment">// Initialize the vector with a function we'd like to interpolate:</span> 205 <span class="keyword">for</span> <span class="special">(</span><span class="identifier">size_t</span> <span class="identifier">i</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="identifier">i</span> <span class="special"><</span> <span class="identifier">v</span><span class="special">.</span><span class="identifier">size</span><span class="special">();</span> <span class="special">++</span><span class="identifier">i</span><span class="special">)</span> 206 <span class="special">{</span> 207 <span class="identifier">v</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sin</span><span class="special">(</span><span class="identifier">i</span><span class="special">*</span><span class="identifier">step</span><span class="special">);</span> 208 <span class="special">}</span> 209 <span class="comment">// We could define an arbitrary start time, but for now we'll just use 0:</span> 210 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">interpolators</span><span class="special">::</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">spline</span><span class="special">(</span><span class="identifier">v</span><span class="special">.</span><span class="identifier">data</span><span class="special">(),</span> <span class="identifier">v</span><span class="special">.</span><span class="identifier">size</span><span class="special">(),</span> <span class="number">0</span> <span class="comment">/* start time */</span><span class="special">,</span> <span class="identifier">step</span><span class="special">);</span> 211 212 <span class="comment">// Now we can evaluate the spline wherever we please.</span> 213 <span class="identifier">std</span><span class="special">::</span><span class="identifier">mt19937</span> <span class="identifier">gen</span><span class="special">;</span> 214 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">random</span><span class="special">::</span><span class="identifier">uniform_real_distribution</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">absissa</span><span class="special">(</span><span class="number">0</span><span class="special">,</span> <span class="identifier">v</span><span class="special">.</span><span class="identifier">size</span><span class="special">()*</span><span class="identifier">step</span><span class="special">);</span> 215 <span class="keyword">for</span> <span class="special">(</span><span class="identifier">size_t</span> <span class="identifier">i</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="identifier">i</span> <span class="special"><</span> <span class="number">10</span><span class="special">;</span> <span class="special">++</span><span class="identifier">i</span><span class="special">)</span> 216 <span class="special">{</span> 217 <span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">absissa</span><span class="special">(</span><span class="identifier">gen</span><span class="special">);</span> 218 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"sin("</span> <span class="special"><<</span> <span class="identifier">x</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">sin</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", spline interpolation gives "</span> <span class="special"><<</span> <span class="identifier">spline</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 219 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cos("</span> <span class="special"><<</span> <span class="identifier">x</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">cos</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", spline derivative interpolation gives "</span> <span class="special"><<</span> <span class="identifier">spline</span><span class="special">.</span><span class="identifier">prime</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 220 <span class="special">}</span> 221 222 <span class="comment">// The next example is less trivial:</span> 223 <span class="comment">// We will try to figure out when the population of the United States crossed 100 million.</span> 224 <span class="comment">// Since the census is taken every 10 years, the data is equally spaced, so we can use the cubic b spline.</span> 225 <span class="comment">// Data taken from https://en.wikipedia.org/wiki/United_States_Census</span> 226 <span class="comment">// We'll start at the year 1860:</span> 227 <span class="keyword">double</span> <span class="identifier">t0</span> <span class="special">=</span> <span class="number">1860</span><span class="special">;</span> 228 <span class="keyword">double</span> <span class="identifier">time_step</span> <span class="special">=</span> <span class="number">10</span><span class="special">;</span> 229 <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">population</span><span class="special">{</span><span class="number">31443321</span><span class="special">,</span> <span class="comment">/* 1860 */</span> 230 <span class="number">39818449</span><span class="special">,</span> <span class="comment">/* 1870 */</span> 231 <span class="number">50189209</span><span class="special">,</span> <span class="comment">/* 1880 */</span> 232 <span class="number">62947714</span><span class="special">,</span> <span class="comment">/* 1890 */</span> 233 <span class="number">76212168</span><span class="special">,</span> <span class="comment">/* 1900 */</span> 234 <span class="number">92228496</span><span class="special">,</span> <span class="comment">/* 1910 */</span> 235 <span class="number">106021537</span><span class="special">,</span> <span class="comment">/* 1920 */</span> 236 <span class="number">122775046</span><span class="special">,</span> <span class="comment">/* 1930 */</span> 237 <span class="number">132164569</span><span class="special">,</span> <span class="comment">/* 1940 */</span> 238 <span class="number">150697361</span><span class="special">,</span> <span class="comment">/* 1950 */</span> 239 <span class="number">179323175</span><span class="special">};/*</span> <span class="number">1960</span> <span class="special">*/</span> 240 241 <span class="comment">// An eyeball estimate indicates that the population crossed 100 million around 1915.</span> 242 <span class="comment">// Let's see what interpolation says:</span> 243 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">interpolators</span><span class="special">::</span><span class="identifier">cardinal_cubic_b_spline</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">p</span><span class="special">(</span><span class="identifier">population</span><span class="special">.</span><span class="identifier">data</span><span class="special">(),</span> <span class="identifier">population</span><span class="special">.</span><span class="identifier">size</span><span class="special">(),</span> <span class="identifier">t0</span><span class="special">,</span> <span class="identifier">time_step</span><span class="special">);</span> 244 245 <span class="comment">// Now create a function which has a zero at p = 100,000,000:</span> 246 <span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[=](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">){</span> <span class="keyword">return</span> <span class="identifier">p</span><span class="special">(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">-</span> <span class="number">100000000</span><span class="special">;</span> <span class="special">};</span> 247 248 <span class="comment">// Boost includes a bisection algorithm, which is robust, though not as fast as some others</span> 249 <span class="comment">// we provide, but let's try that first. We need a termination condition for it, which</span> 250 <span class="comment">// takes the two endpoints of the range and returns either true (stop) or false (keep going),</span> 251 <span class="comment">// we could use a predefined one such as boost::math::tools::eps_tolerance<double>, but that</span> 252 <span class="comment">// won't stop until we have full double precision which is overkill, since we just need the</span> 253 <span class="comment">// endpoint to yield the same month. While we're at it, we'll keep track of the number of</span> 254 <span class="comment">// iterations required too, though this is strictly optional:</span> 255 256 <span class="keyword">auto</span> <span class="identifier">termination</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">left</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">right</span><span class="special">)</span> 257 <span class="special">{</span> 258 <span class="keyword">double</span> <span class="identifier">left_month</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">round</span><span class="special">((</span><span class="identifier">left</span> <span class="special">-</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">left</span><span class="special">))</span> <span class="special">*</span> <span class="number">12</span> <span class="special">+</span> <span class="number">1</span><span class="special">);</span> 259 <span class="keyword">double</span> <span class="identifier">right_month</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">round</span><span class="special">((</span><span class="identifier">right</span> <span class="special">-</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">right</span><span class="special">))</span> <span class="special">*</span> <span class="number">12</span> <span class="special">+</span> <span class="number">1</span><span class="special">);</span> 260 <span class="keyword">return</span> <span class="special">(</span><span class="identifier">left_month</span> <span class="special">==</span> <span class="identifier">right_month</span><span class="special">)</span> <span class="special">&&</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">left</span><span class="special">)</span> <span class="special">==</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">right</span><span class="special">));</span> 261 <span class="special">};</span> 262 <span class="identifier">std</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">iterations</span> <span class="special">=</span> <span class="number">1000</span><span class="special">;</span> 263 <span class="keyword">auto</span> <span class="identifier">result</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">bisect</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">1910.0</span><span class="special">,</span> <span class="number">1920.0</span><span class="special">,</span> <span class="identifier">termination</span><span class="special">,</span> <span class="identifier">iterations</span><span class="special">);</span> 264 <span class="keyword">auto</span> <span class="identifier">time</span> <span class="special">=</span> <span class="identifier">result</span><span class="special">.</span><span class="identifier">first</span><span class="special">;</span> <span class="comment">// termination condition ensures that both endpoints yield the same result</span> 265 <span class="keyword">auto</span> <span class="identifier">month</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">round</span><span class="special">((</span><span class="identifier">time</span> <span class="special">-</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">time</span><span class="special">))*</span><span class="number">12</span> <span class="special">+</span> <span class="number">1</span><span class="special">);</span> 266 <span class="keyword">auto</span> <span class="identifier">year</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">time</span><span class="special">);</span> 267 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The population of the United States surpassed 100 million on the "</span><span class="special">;</span> 268 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">month</span> <span class="special"><<</span> <span class="string">"th month of "</span> <span class="special"><<</span> <span class="identifier">year</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 269 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Found in "</span> <span class="special"><<</span> <span class="identifier">iterations</span> <span class="special"><<</span> <span class="string">" iterations"</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 270 271 <span class="comment">// Since the cubic B spline offers the first derivative, we could equally have used Newton iterations,</span> 272 <span class="comment">// this takes "number of bits correct" as a termination condition - 20 should be plenty for what we need,</span> 273 <span class="comment">// and once again, we track how many iterations are taken:</span> 274 275 <span class="keyword">auto</span> <span class="identifier">f_n</span> <span class="special">=</span> <span class="special">[=](</span><span class="keyword">double</span> <span class="identifier">t</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span><span class="identifier">p</span><span class="special">(</span><span class="identifier">t</span><span class="special">)</span> <span class="special">-</span> <span class="number">100000000</span><span class="special">,</span> <span class="identifier">p</span><span class="special">.</span><span class="identifier">prime</span><span class="special">(</span><span class="identifier">t</span><span class="special">));</span> <span class="special">};</span> 276 <span class="identifier">iterations</span> <span class="special">=</span> <span class="number">1000</span><span class="special">;</span> 277 <span class="identifier">time</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">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">f_n</span><span class="special">,</span> <span class="number">1910.0</span><span class="special">,</span> <span class="number">1900.0</span><span class="special">,</span> <span class="number">2000.0</span><span class="special">,</span> <span class="number">20</span><span class="special">,</span> <span class="identifier">iterations</span><span class="special">);</span> 278 <span class="identifier">month</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">round</span><span class="special">((</span><span class="identifier">time</span> <span class="special">-</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">time</span><span class="special">))*</span><span class="number">12</span> <span class="special">+</span> <span class="number">1</span><span class="special">);</span> 279 <span class="identifier">year</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">floor</span><span class="special">(</span><span class="identifier">time</span><span class="special">);</span> 280 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"The population of the United States surpassed 100 million on the "</span><span class="special">;</span> 281 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">month</span> <span class="special"><<</span> <span class="string">"th month of "</span> <span class="special"><<</span> <span class="identifier">year</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 282 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Found in "</span> <span class="special"><<</span> <span class="identifier">iterations</span> <span class="special"><<</span> <span class="string">" iterations"</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 283 284<span class="special">}</span> 285</pre> 286<pre class="programlisting"><span class="identifier">Program</span> <span class="identifier">output</span> <span class="identifier">is</span><span class="special">:</span> 287</pre> 288<pre class="programlisting">sin(4.07362) = -0.802829, spline interpolation gives - 0.802829 289cos(4.07362) = -0.596209, spline derivative interpolation gives - 0.596209 290sin(0.677385) = 0.626758, spline interpolation gives 0.626758 291cos(0.677385) = 0.779214, spline derivative interpolation gives 0.779214 292sin(4.52896) = -0.983224, spline interpolation gives - 0.983224 293cos(4.52896) = -0.182402, spline derivative interpolation gives - 0.182402 294sin(4.17504) = -0.85907, spline interpolation gives - 0.85907 295cos(4.17504) = -0.511858, spline derivative interpolation gives - 0.511858 296sin(0.634934) = 0.593124, spline interpolation gives 0.593124 297cos(0.634934) = 0.805111, spline derivative interpolation gives 0.805111 298sin(4.84434) = -0.991307, spline interpolation gives - 0.991307 299cos(4.84434) = 0.131567, spline derivative interpolation gives 0.131567 300sin(4.56688) = -0.989432, spline interpolation gives - 0.989432 301cos(4.56688) = -0.144997, spline derivative interpolation gives - 0.144997 302sin(1.10517) = 0.893541, spline interpolation gives 0.893541 303cos(1.10517) = 0.448982, spline derivative interpolation gives 0.448982 304sin(3.1618) = -0.0202022, spline interpolation gives - 0.0202022 305cos(3.1618) = -0.999796, spline derivative interpolation gives - 0.999796 306sin(1.54084) = 0.999551, spline interpolation gives 0.999551 307cos(1.54084) = 0.0299566, spline derivative interpolation gives 0.0299566 308The population of the United States surpassed 100 million on the 11th month of 1915 309Found in 12 iterations 310The population of the United States surpassed 100 million on the 11th month of 1915 311Found in 3 iterations 312</pre> 313</div> 314<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 315<td align="left"></td> 316<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 317 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 318 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 319 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 320 Daryle Walker and Xiaogang Zhang<p> 321 Distributed under the Boost Software License, Version 1.0. (See accompanying 322 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>) 323 </p> 324</div></td> 325</tr></table> 326<hr> 327<div class="spirit-nav"> 328<a accesskey="p" href="../interpolation.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../interpolation.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="cardinal_quadratic_b.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 329</div> 330</body> 331</html> 332
