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₀ 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₁ := B₀ ∗ B₀, and B_n+1 := B_n ∗ B₀. For example, B₁(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);

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₁ is what Schoenberg calls M₂. 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].

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.

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.

Accuracy

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

References