diff --git a/doc/graphs/central_b_spline_derivatives.svg b/doc/graphs/central_b_spline_derivatives.svg
new file mode 100644
index 000000000..490f93fd5
--- /dev/null
+++ b/doc/graphs/central_b_spline_derivatives.svg
@@ -0,0 +1,56 @@
+
+
diff --git a/doc/graphs/central_b_spline_second_derivatives.svg b/doc/graphs/central_b_spline_second_derivatives.svg
new file mode 100644
index 000000000..7bd209f04
--- /dev/null
+++ b/doc/graphs/central_b_spline_second_derivatives.svg
@@ -0,0 +1,54 @@
+
+
diff --git a/doc/sf/cardinal_b_splines.qbk b/doc/sf/cardinal_b_splines.qbk
index 7f7c4e15e..79e61c49a 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,14 +40,22 @@ 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]
+[$../graphs/central_b_spline_derivatives.svg]
+[$../graphs/central_b_spline_second_derivatives.svg]
[h3 Caveats]
diff --git a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
index 1347cbc30..684fb76b9 100644
--- a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
+++ b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
@@ -91,7 +91,7 @@ Real b3_spline_double_prime(Real x)
if (x < 1)
{
- return (3*boost::math::constants::half()*x - 2) + x*(3*boost::math::constants::half());
+ return 3*x - 2;
}
if (x < 2)
{
diff --git a/include/boost/math/special_functions/cardinal_b_spline.hpp b/include/boost/math/special_functions/cardinal_b_spline.hpp
index 399934ffd..dfb4cf835 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,185 @@ 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);
+ }
+ }
+
+ if (n==1)
+ {
+ if (x==0)
+ {
+ return Real(0);
+ }
+ if (x==1)
+ {
+ return -Real(1)/Real(2);
+ }
+ return Real(-1);
+ }
+
+
+ Real supp_max = (n+1)/Real(2);
+ if (x >= supp_max)
+ {
+ return Real(0);
+ }
+
+ // 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 cardinal_b_spline_double_prime(Real x)
+{
+ static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integer types.");
+ static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required.");
+
+ if (x < 0)
+ {
+ // All B-splines are even functions, so second derivatives are even:
+ return cardinal_b_spline_double_prime(-x);
+ }
+
+
+ Real supp_max = (n+1)/Real(2);
+ if (x >= supp_max)
+ {
+ return Real(0);
+ }
+
+ // 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 - 2; ++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[2] - 2*v[1] + v[0];
+}
+
+
template
-Real forward_cardinal_b_spline(Real x) {
+Real forward_cardinal_b_spline(Real x)
+{
+ static_assert(!std::is_integral::value, "Cardinal B-splines do not work with integral types.");
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..0ec92cc95 100644
--- a/test/cardinal_b_spline_test.cpp
+++ b/test/cardinal_b_spline_test.cpp
@@ -10,8 +10,8 @@
#include
#include
#include
-#include
#include
+#include
#ifdef BOOST_HAS_FLOAT128
#include
using boost::multiprecision::float128;
@@ -19,7 +19,10 @@ 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;
+using boost::math::cardinal_b_spline_double_prime;
+
template
void test_box()
@@ -27,6 +30,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 +41,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 +71,23 @@ void test_hat()
{
std::cerr << " Problem at t = " << t << "\n";
}
+ if (t == -Real(1)) {
+ if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0)) {
+ std::cout << " Problem at t = " << t << "\n";
+ }
+ }
+ else if (t == Real(1)) {
+ CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0);
+ }
+ else 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)
@@ -92,10 +115,33 @@ void test_quadratic()
return (2-t1*t1 -t2*t2)/2;
};
+ auto b2_prime = [&](Real x)->Real {
+ Real absx = abs(x);
+ Real signx = 1;
+ if (x < 0) {
+ signx = -1;
+ }
+ if (absx >= 3/Real(2)) {
+ return Real(0);
+ }
+ if (absx >= 1/Real(2)) {
+ return (absx - 3/Real(2))*signx;
+ }
+ return -2*absx*signx;
+ };
+
+
Real h = 1/Real(256);
for (Real t = -5; t <= 5; t += h) {
Real expected = b2(t);
CHECK_ULP_CLOSE(expected, cardinal_b_spline<2>(t), 0);
+ expected = b2_prime(t);
+
+ if (!CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<2>(t), 0))
+ {
+ std::cerr << " Problem at t = " << t << "\n";
+ }
+
}
}
@@ -113,6 +159,18 @@ void test_cubic()
expected = Real(0);
computed = cardinal_b_spline<3, Real>(2);
CHECK_ULP_CLOSE(expected, computed, 0);
+
+ Real h = 1/Real(256);
+ for (Real t = -4; t <= 4; t += h) {
+ expected = boost::math::detail::b3_spline_prime(t);
+ computed = cardinal_b_spline_prime<3>(t);
+ CHECK_ULP_CLOSE(expected, computed, 0);
+ expected = boost::math::detail::b3_spline_double_prime(t);
+ computed = cardinal_b_spline_double_prime<3>(t);
+ if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
+ std::cerr << " Problem at t = " << t << "\n";
+ }
+ }
}
template
@@ -136,6 +194,23 @@ void test_quintic()
}
+template
+void test_b_spline_derivatives()
+{
+ Real h = 1/Real(256);
+ Real supp = (n+Real(1))/Real(2);
+ for (Real t = -supp - 1; t <= supp+1; t+= h)
+ {
+ Real expected = cardinal_b_spline(t+Real(1)/Real(2)) - cardinal_b_spline(t - Real(1)/Real(2));
+ Real computed = cardinal_b_spline_prime(t);
+ CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits::epsilon());
+
+ expected = cardinal_b_spline(t+1) - 2*cardinal_b_spline(t) + cardinal_b_spline(t-1);
+ computed = cardinal_b_spline_double_prime(t);
+ CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits::epsilon());
+ }
+}
+
template
void test_partition_of_unity()
{
@@ -186,6 +261,23 @@ int main()
test_partition_of_unity<5, double>();
test_partition_of_unity<6, double>();
+ test_b_spline_derivatives<3, double>();
+ test_b_spline_derivatives<4, double>();
+ test_b_spline_derivatives<5, double>();
+ test_b_spline_derivatives<6, double>();
+ test_b_spline_derivatives<7, double>();
+ test_b_spline_derivatives<8, double>();
+ test_b_spline_derivatives<9, double>();
+
+ test_b_spline_derivatives<3, long double>();
+ test_b_spline_derivatives<4, long double>();
+ test_b_spline_derivatives<5, long double>();
+ test_b_spline_derivatives<6, long double>();
+ test_b_spline_derivatives<7, long double>();
+ test_b_spline_derivatives<8, long double>();
+ test_b_spline_derivatives<9, long double>();
+
+
#ifdef BOOST_HAS_FLOAT128
test_box();
test_hat();