From e327be5887a5fd8aa4ae3544b6738bdb50d43f10 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 11 Aug 2019 14:52:47 -0400 Subject: [PATCH] Cardinal B-spline derivatives [CI SKIP] --- doc/sf/cardinal_b_splines.qbk | 14 +- .../special_functions/cardinal_b_spline.hpp | 172 +++++++++++++----- test/cardinal_b_spline_test.cpp | 13 ++ 3 files changed, 148 insertions(+), 51 deletions(-) diff --git a/doc/sf/cardinal_b_splines.qbk b/doc/sf/cardinal_b_splines.qbk index 7f7c4e15e..f4e2f3938 100644 --- a/doc/sf/cardinal_b_splines.qbk +++ b/doc/sf/cardinal_b_splines.qbk @@ -16,7 +16,13 @@ namespace boost{ namespace math{ template - Real cardinal_b_spline(Real x); + auto cardinal_b_spline(Real x); + + template + auto cardinal_b_spline_prime(Real x); + + template + auto cardinal_b_spline_double_prime(Real x); template Real forward_cardinal_b_spline(Real x); @@ -34,11 +40,17 @@ For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere 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] diff --git a/include/boost/math/special_functions/cardinal_b_spline.hpp b/include/boost/math/special_functions/cardinal_b_spline.hpp index 399934ffd..4e40058ed 100644 --- a/include/boost/math/special_functions/cardinal_b_spline.hpp +++ b/include/boost/math/special_functions/cardinal_b_spline.hpp @@ -8,6 +8,8 @@ #include #include +#include +#include namespace boost { namespace math { @@ -30,65 +32,135 @@ namespace detail { template Real cardinal_b_spline(Real x) { - if (x < 0) { - // All B-splines are even functions: - return cardinal_b_spline(-x); - } + static_assert(!std::is_integral::value, "Does not work with integral types."); - // Someday we'll be able to 'if constexpr' on many compilers. - // Not today though! In any case, it's not a big performance hit. - if /*constexpr*/ (n==0) - { - if (x < Real(1)/Real(2)) { - return Real(1); + if (x < 0) { + // All B-splines are even functions: + return cardinal_b_spline(-x); } - else if (x == Real(1)/Real(2)) { - return Real(1)/Real(2); - } - else { - return Real(0); - } - } - if /*constexpr*/ (n==1) { - return detail::B1(x); - } - - Real supp_max = (n+1)/Real(2); - if (x >= supp_max) - { - return Real(0); - } - - // Fill v with values of B1: - // At most two of these terms are nonzero, and at least 1. - // There is only one non-zero term when n is odd and x = 0. - std::array v; - Real z = x + 1 - supp_max; - for (unsigned i = 0; i < n; ++i) - { - v[i] = detail::B1(z); - z += 1; - } - - Real smx = supp_max - x; - for (unsigned j = 2; j <= n; ++j) - { - Real a = (j + 1 - smx); - Real b = smx; - for(unsigned k = 0; k <= n - j; ++k) + if (n==0) { - v[k] = (a*v[k+1] + b*v[k])/Real(j); - a += 1; - b -= 1; + if (x < Real(1)/Real(2)) { + return Real(1); + } + else if (x == Real(1)/Real(2)) { + return Real(1)/Real(2); + } + else { + return Real(0); + } } - } - return v[0]; + if (n==1) + { + return detail::B1(x); + } + + Real supp_max = (n+1)/Real(2); + if (x >= supp_max) + { + return Real(0); + } + + // Fill v with values of B1: + // At most two of these terms are nonzero, and at least 1. + // There is only one non-zero term when n is odd and x = 0. + std::array v; + Real z = x + 1 - supp_max; + for (unsigned i = 0; i < n; ++i) + { + v[i] = detail::B1(z); + z += 1; + } + + Real smx = supp_max - x; + for (unsigned j = 2; j <= n; ++j) + { + Real a = (j + 1 - smx); + Real b = smx; + for(unsigned k = 0; k <= n - j; ++k) + { + v[k] = (a*v[k+1] + b*v[k])/Real(j); + a += 1; + b -= 1; + } + } + + return v[0]; +} + + +template +Real cardinal_b_spline_prime(Real x) +{ + static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integer types."); + + if (x < 0) + { + // All B-splines are even functions, so derivatives are odd: + return -cardinal_b_spline_prime(-x); + } + + + if (n==0) + { + // Kinda crazy but you get what you ask for! + if (x == Real(1)/Real(2)) + { + return std::numeric_limits::infinity(); + } + else + { + return Real(0); + } + } + + Real supp_max = (n+1)/Real(2); + if (x >= supp_max) + { + return Real(0); + } + + if (n==1) { + if (x==0) { + return Real(0); + } + return Real(-1); + } + + // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2): + std::array v; + Real z = x + 1 - supp_max; + for (unsigned i = 0; i < n; ++i) + { + v[i] = detail::B1(z); + z += 1; + } + + Real smx = supp_max - x; + for (unsigned j = 2; j <= n - 1; ++j) + { + Real a = (j + 1 - smx); + Real b = smx; + for(unsigned k = 0; k <= n - j; ++k) + { + v[k] = (a*v[k+1] + b*v[k])/Real(j); + a += 1; + b -= 1; + } + } + + return v[1] - v[0]; } template -Real forward_cardinal_b_spline(Real x) { +auto forward_cardinal_b_spline(Real x) +{ + if constexpr (std::is_integral::value) + { + return forward_cardinal_b_spline(x); + } return cardinal_b_spline(x - (n+1)/Real(2)); } diff --git a/test/cardinal_b_spline_test.cpp b/test/cardinal_b_spline_test.cpp index b819abba6..20c5c5dc7 100644 --- a/test/cardinal_b_spline_test.cpp +++ b/test/cardinal_b_spline_test.cpp @@ -19,6 +19,7 @@ using boost::multiprecision::float128; using std::abs; using boost::math::cardinal_b_spline; +using boost::math::cardinal_b_spline_prime; using boost::math::forward_cardinal_b_spline; template @@ -27,6 +28,7 @@ void test_box() Real t = cardinal_b_spline<0>(Real(1.1)); Real expected = 0; CHECK_ULP_CLOSE(expected, t, 0); + CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0); t = cardinal_b_spline<0>(Real(-1.1)); expected = 0; @@ -37,6 +39,8 @@ void test_box() { expected = 1; CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0); + expected = 0; + CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0); } for (t = h; t < 1; t += h) @@ -65,6 +69,15 @@ void test_hat() { std::cerr << " Problem at t = " << t << "\n"; } + if (t < 0) { + CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0); + } + else if (t == 0) { + CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime<1>(t), 0); + } + else if (t > 0) { + CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime<1>(t), 0); + } } for (t = 0; t < 2; t += h)