1[/ 2 Copyright 2019, Nick Thompson 3 Distributed under the Boost Software License, Version 1.0. 4 (See accompanying file LICENSE_1_0.txt or copy at 5 http://www.boost.org/LICENSE_1_0.txt). 6] 7 8[section:cardinal_b_splines Cardinal B-splines] 9 10[h4 Synopsis] 11 12`` 13#include <boost/math/special_functions/cardinal_b_spline.hpp> 14`` 15 16 namespace boost{ namespace math{ 17 18 template<unsigned n, typename Real> 19 auto cardinal_b_spline(Real x); 20 21 template<unsigned n, typename Real> 22 auto cardinal_b_spline_prime(Real x); 23 24 template<unsigned n, typename Real> 25 auto cardinal_b_spline_double_prime(Real x); 26 27 template<unsigned n, typename Real> 28 Real forward_cardinal_b_spline(Real x); 29 30 }} // namespaces 31 32Cardinal /B/-splines are a family of compactly supported functions useful for the smooth interpolation of tables. 33 34The first /B/-spline /B/[sub 0] is simply a box function: 35It takes the value one inside the interval \[-1/2, 1/2\], and is zero elsewhere. 36/B/-splines of higher smoothness are constructed by iterative convolution, namely, /B/[sub 1] := /B/[sub 0] \u2217 /B/[sub 0], 37and /B/[sub /n/+1] := /B/[sub /n/ ] \u2217 /B/[sub 0]. 38For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere, so it is a hat function. 39 40A basic usage is as follows: 41 42 using boost::math::cardinal_b_spline; 43 using boost::math::cardinal_b_spline_prime; 44 using boost::math::cardinal_b_spline_double_prime; 45 double x = 0.5; 46 // B₀(x), the box function: 47 double y = cardinal_b_spline<0>(x); 48 // B₁(x), the hat function: 49 y = cardinal_b_spline<1>(x); 50 // First derivative of B₃: 51 yp = cardinal_b_spline_prime<3>(x); 52 // Second derivative of B₃: 53 ypp = cardinal_b_spline_double_prime<3>(x); 54 55 56[$../graphs/central_b_splines.svg] 57[$../graphs/central_b_spline_derivatives.svg] 58[$../graphs/central_b_spline_second_derivatives.svg] 59 60[h3 Caveats] 61 62Numerous notational conventions for /B/-splines exist. 63Whereas Boost.Math (following Kress) zero indexes the /B/-splines, 64other authors (such as Schoenberg and Massopust) use 1-based indexing. 65So (for example) Boost.Math's hat function /B/[sub 1] is what Schoenberg calls /M/[sub 2]. 66Mathematica, like Boost, uses the zero-indexing convention for its `BSplineCurve`. 67 68Even the support of the splines is not agreed upon. 69Mathematica starts the support of the splines at zero and rescales the independent variable so that the support of every member is \[0, 1\]. 70Massopust as well as Unser puts the support of the /B/-splines at \[0, /n/\], whereas Kress centers them at zero. 71Schoenberg distinguishes between the two cases, called the splines starting at zero forward splines, and the ones symmetric about zero /central/. 72 73The /B/-splines of Boost.Math are central, with support support \[-(/n/+1)\/2, (/n/+1)\/2\]. 74If necessary, the forward splines can be evaluated by using `forward_cardinal_b_spline`, whose support is \[0, /n/+1\]. 75 76[h3 Implementation] 77 78The 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. 79Benchmarks revealed that the time to calculate the indexes of the rhombohedral array exceed the time to simply add zeroes together, except for /n/ > 18. 80Since few people use /B/ splines of degree 18, the triangular array is used. 81 82[h3 Performance] 83 84Double precision timing on a consumer x86 laptop is shown below: 85 86`` 87Run on (16 X 4300 MHz CPU s) 88CPU Caches: 89 L1 Data 32K (x8) 90 L1 Instruction 32K (x8) 91 L2 Unified 1024K (x8) 92 L3 Unified 11264K (x1) 93Load Average: 0.21, 0.33, 0.29 94----------------------------------------- 95Benchmark Time 96----------------------------------------- 97CardinalBSpline<1, double> 18.8 ns 98CardinalBSpline<2, double> 25.3 ns 99CardinalBSpline<3, double> 29.3 ns 100CardinalBSpline<4, double> 33.8 ns 101CardinalBSpline<5, double> 36.7 ns 102CardinalBSpline<6, double> 39.1 ns 103CardinalBSpline<7, double> 43.6 ns 104CardinalBSpline<8, double> 62.8 ns 105CardinalBSpline<9, double> 70.2 ns 106CardinalBSpline<10, double> 83.8 ns 107CardinalBSpline<11, double> 94.3 ns 108CardinalBSpline<12, double> 108 ns 109CardinalBSpline<13, double> 122 ns 110CardinalBSpline<14, double> 138 ns 111CardinalBSpline<15, double> 155 ns 112CardinalBSpline<16, double> 170 ns 113CardinalBSpline<17, double> 192 ns 114CardinalBSpline<18, double> 174 ns 115CardinalBSpline<19, double> 180 ns 116CardinalBSpline<20, double> 194 ns 117UniformReal<double> 11.5 ns 118`` 119 120A 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. 121 122[h3 Accuracy] 123 124Some representative ULP plots are shown below. 125The error grows linearly with /n/, as expected from Cox equation 10.5. 126 127[$../graphs/b_spline_ulp_3.svg] 128[$../graphs/b_spline_ulp_5.svg] 129[$../graphs/b_spline_ulp_9.svg] 130 131[h3 References] 132 133* I.J. Schoenberg, ['Cardinal Spline Interpolation], SIAM Volume 12, 1973 134* Rainer Kress, ['Numerical Analysis], Springer, 1998 135* Peter Massopust, ['On Some Generalizations of B-splines], arxiv preprint, 2019 136* Michael Unser and Thierry Blu, ['Fractional Splines and Wavelets], SIAM Review 2000, Volume 42, No. 1 137* Cox, Maurice G. ['The numerical evaluation of B-splines.], IMA Journal of Applied Mathematics 10.2 (1972): 134-149. 138 139[endsect] 140
