EQ2EMu/EQ2/source/depends/boost_1_72_0/libs/math/doc/sf/cardinal_b_splines.qbk

[/
  Copyright 2019, Nick Thompson
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]

[section:cardinal_b_splines Cardinal B-splines]

[h4 Synopsis]

``
#include <boost/math/special_functions/cardinal_b_spline.hpp>
``

   namespace boost{ namespace math{

   template<unsigned n, typename Real>
   auto cardinal_b_spline(Real x);

   template<unsigned n, typename Real>
   auto cardinal_b_spline_prime(Real x);

   template<unsigned n, typename Real>
   auto cardinal_b_spline_double_prime(Real x);

   template<unsigned n, typename Real>
   Real forward_cardinal_b_spline(Real x);

   }} // namespaces

Cardinal /B/-splines are a family of compactly supported functions useful for the smooth interpolation of tables.

The first /B/-spline /B/[sub 0] is simply a box function:
It takes the value one inside the interval \[-1/2, 1/2\], and is zero elsewhere.
/B/-splines of higher smoothness are constructed by iterative convolution, namely, /B/[sub 1] := /B/[sub 0] \u2217 /B/[sub 0],
and /B/[sub /n/+1] := /B/[sub /n/ ] \u2217 /B/[sub 0].
For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere, so it is a hat function.

A basic usage is as follows:

    using boost::math::cardinal_b_spline;
    using boost::math::cardinal_b_spline_prime;
    using boost::math::cardinal_b_spline_double_prime;
    double x = 0.5;
    // B₀(x), the box function:
    double y = cardinal_b_spline<0>(x);
    // B₁(x), the hat function:
    y = cardinal_b_spline<1>(x);
    // First derivative of B₃:
    yp = cardinal_b_spline_prime<3>(x);
    // Second derivative of B₃:
    ypp = cardinal_b_spline_double_prime<3>(x);


[$../graphs/central_b_splines.svg]
[$../graphs/central_b_spline_derivatives.svg]
[$../graphs/central_b_spline_second_derivatives.svg]

[h3 Caveats]

Numerous notational conventions for /B/-splines exist.
Whereas Boost.Math (following Kress) zero indexes the /B/-splines,
other authors (such as Schoenberg and Massopust) use 1-based indexing.
So (for example) Boost.Math's hat function /B/[sub 1] is what Schoenberg calls /M/[sub 2].
Mathematica, like Boost, uses the zero-indexing convention for its `BSplineCurve`.

Even the support of the splines is not agreed upon.
Mathematica starts the support of the splines at zero and rescales the independent variable so that the support of every member is \[0, 1\].
Massopust as well as Unser puts the support of the /B/-splines at \[0, /n/\], whereas Kress centers them at zero.
Schoenberg distinguishes between the two cases, called the splines starting at zero forward splines, and the ones symmetric about zero /central/.

The /B/-splines of Boost.Math are central, with support support \[-(/n/+1)\/2, (/n/+1)\/2\].
If necessary, the forward splines can be evaluated by using `forward_cardinal_b_spline`, whose support is \[0, /n/+1\].

[h3 Implementation]

The implementation follows Maurice Cox' 1972 paper 'The Numerical Evaluation of B-splines', and uses the triangular array of Algorithm 6.1 of the reference rather than the rhombohedral array of Algorithm 6.2.
Benchmarks revealed that the time to calculate the indexes of the rhombohedral array exceed the time to simply add zeroes together, except for /n/ > 18.
Since few people use /B/ splines of degree 18, the triangular array is used.

[h3 Performance]

Double precision timing on a consumer x86 laptop is shown below:

``
Run on (16 X 4300 MHz CPU s)
CPU Caches:
  L1 Data 32K (x8)
  L1 Instruction 32K (x8)
  L2 Unified 1024K (x8)
  L3 Unified 11264K (x1)
Load Average: 0.21, 0.33, 0.29
-----------------------------------------
Benchmark                            Time
-----------------------------------------
CardinalBSpline<1, double>        18.8 ns
CardinalBSpline<2, double>        25.3 ns
CardinalBSpline<3, double>        29.3 ns
CardinalBSpline<4, double>        33.8 ns
CardinalBSpline<5, double>        36.7 ns
CardinalBSpline<6, double>        39.1 ns
CardinalBSpline<7, double>        43.6 ns
CardinalBSpline<8, double>        62.8 ns
CardinalBSpline<9, double>        70.2 ns
CardinalBSpline<10, double>       83.8 ns
CardinalBSpline<11, double>       94.3 ns
CardinalBSpline<12, double>        108 ns
CardinalBSpline<13, double>        122 ns
CardinalBSpline<14, double>        138 ns
CardinalBSpline<15, double>        155 ns
CardinalBSpline<16, double>        170 ns
CardinalBSpline<17, double>        192 ns
CardinalBSpline<18, double>        174 ns
CardinalBSpline<19, double>        180 ns
CardinalBSpline<20, double>        194 ns
UniformReal<double>               11.5 ns
``

A uniformly distributed random number within the support of the spline is generated for the argument, so subtracting 11.5 ns from these gives a good idea of the performance of the calls.

[h3 Accuracy]

Some representative ULP plots are shown below.
The error grows linearly with /n/, as expected from Cox equation 10.5.

[$../graphs/b_spline_ulp_3.svg]
[$../graphs/b_spline_ulp_5.svg]
[$../graphs/b_spline_ulp_9.svg]

[h3 References]

* I.J. Schoenberg, ['Cardinal Spline Interpolation], SIAM Volume 12, 1973
* Rainer Kress, ['Numerical Analysis], Springer, 1998
* Peter Massopust, ['On Some Generalizations of B-splines], arxiv preprint, 2019
* Michael Unser and Thierry Blu, ['Fractional Splines and Wavelets], SIAM Review 2000, Volume 42, No. 1
* Cox, Maurice G. ['The numerical evaluation of B-splines.], IMA Journal of Applied Mathematics 10.2 (1972): 134-149.

[endsect]