boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-25 16:32:15 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #240 from boostorg/cardinal_b_spline_derivatives

Browse Source 
Cardinal B-spline derivatives
This commit is contained in:
Nick
2019-08-18 09:21:41 -04:00
committed by GitHub
parent 17faf9ba06 3291ac891b
commit cf48e191a6
6 changed files with 391 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
test/cardinal_b_spline_test.cpp
View File

@@ -10,8 +10,8 @@
#include <utility>
#include <random>
#include <cmath>
#include <boost/core/demangle.hpp>
#include <boost/math/special_functions/cardinal_b_spline.hpp>
#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
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<class Real>
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<Real>(t);
computed = cardinal_b_spline_prime<3>(t);
CHECK_ULP_CLOSE(expected, computed, 0);
expected = boost::math::detail::b3_spline_double_prime<Real>(t);
computed = cardinal_b_spline_double_prime<3>(t);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Problem at t = " << t << "\n";
}
}
}
template<class Real>
@@ -136,6 +194,23 @@ void test_quintic()
}
template<unsigned n, typename Real>
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<n-1>(t+Real(1)/Real(2)) - cardinal_b_spline<n-1>(t - Real(1)/Real(2));
Real computed = cardinal_b_spline_prime<n>(t);
CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon());
expected = cardinal_b_spline<n-2>(t+1) - 2*cardinal_b_spline<n-2>(t) + cardinal_b_spline<n-2>(t-1);
computed = cardinal_b_spline_double_prime<n>(t);
CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon());
}
}
template<unsigned n, typename Real>
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<float128>();
test_hat<float128>();