libs/math/doc/html/math_toolkit/cubic_b.html

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