boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix septic hermite evaluation.

Browse Source
This commit is contained in:
Nick
2020-03-22 13:03:48 -04:00
parent 2af33f2725
commit fb26bc0cf1
3 changed files with 101 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
include/boost/math/interpolators/detail/septic_hermite_detail.hpp
View File

@@ -184,6 +184,11 @@ public:
return 5*x_.size()*sizeof(Real) + 5*sizeof(x_);
}
std::pair<Real, Real> domain() const
{
return {x_.front(), x_.back()};
}
private:
RandomAccessContainer x_;
RandomAccessContainer y_;
@@ -260,6 +265,9 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s3 - ii;
if (t == 0) {
return y_[i];
}
// See:
// http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
Real t2 = t*t;
@@ -314,6 +322,10 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s3 - ii;
if (t==0)
{
return dy_[i]/inv_dx_;
}
Real y0 = y_[i];
Real y1 = y_[i+1];
@@ -366,6 +378,10 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s3 - ii;
if (t==0)
{
return d2y_[i]*2*inv_dx_*inv_dx_;
}
Real y0 = y_[i];
Real y1 = y_[i+1];
@@ -398,6 +414,11 @@ public:
return 4*y_.size()*sizeof(Real) + 2*sizeof(Real) + 4*sizeof(y_);
}
std::pair<Real, Real> domain() const
{
return {x0_, x0_ + (y_.size()-1)/inv_dx_};
}
private:
RandomAccessContainer y_;
RandomAccessContainer dy_;
@@ -456,6 +477,10 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s3 - ii;
if (t==0)
{
return data_[i][0];
}
Real t2 = t*t;
Real t3 = t2*t;
Real t4 = t3*t;
@@ -508,6 +533,10 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s3 - ii;
if (t == 0)
{
return data_[i][1]*inv_dx_;
}
Real y0 = data_[i][0];
Real y1 = data_[i+1][0];
@@ -560,7 +589,10 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s3 - ii;
if (t == 0)
{
return data_[i][2]*2*inv_dx_*inv_dx_;
}
Real y0 = data_[i][0];
Real y1 = data_[i+1][0];
Real dy0 = data_[i][1];
@@ -592,6 +624,11 @@ public:
return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real) + sizeof(data_);
}
std::pair<Real, Real> domain() const
{
return {x0_, x0_ + (data_.size() -1)/inv_dx_};
}
private:
RandomAccessContainer data_;
Real x0_;

15
include/boost/math/interpolators/septic_hermite.hpp
View File

@@ -44,6 +44,11 @@ public:
return impl_->bytes() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::septic_hermite_detail<RandomAccessContainer>> impl_;
};
@@ -79,6 +84,11 @@ public:
return impl_->bytes() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::cardinal_septic_hermite_detail<RandomAccessContainer>> impl_;
};
@@ -113,6 +123,11 @@ public:
return impl_.size() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::cardinal_septic_hermite_detail_aos<RandomAccessContainer>> impl_;
};

48
test/septic_hermite_test.cpp
View File

@@ -12,6 +12,7 @@
#include <boost/random/uniform_real.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/math/interpolators/septic_hermite.hpp>
#include <boost/math/special_functions/next.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
@@ -74,6 +75,29 @@ void test_constant()
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(t), 24);
}
// Now check the boundaries:
auto [tlo, thi] = csh.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(7), csh(tlo), 2);
CHECK_ULP_CLOSE(Real(7), csh(thi), 2);
CHECK_ULP_CLOSE(Real(7), csh_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(7), csh_aos(thi), 2);
CHECK_ULP_CLOSE(Real(0), csh.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), csh.prime(thi), 2);
CHECK_ULP_CLOSE(Real(0), csh_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), csh_aos.prime(thi), 2);
CHECK_ULP_CLOSE(Real(0), csh.double_prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), csh.double_prime(thi), 2);
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(thi), 2);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
@@ -151,6 +175,30 @@ void test_linear()
CHECK_ULP_CLOSE(Real(1), csh_aos.prime(t), 15);
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(t), 15);
}
// Now check the boundaries:
auto [tlo, thi] = csh.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(tlo), csh(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), csh(thi), 2);
CHECK_ULP_CLOSE(Real(tlo), csh_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), csh_aos(thi), 2);
CHECK_ULP_CLOSE(Real(1), csh.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), csh.prime(thi), 500);
CHECK_ULP_CLOSE(Real(1), csh_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), csh_aos.prime(thi), 500);
CHECK_MOLLIFIED_CLOSE(Real(0), csh.double_prime(tlo), std::numeric_limits<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), csh.double_prime(thi), 400*std::numeric_limits<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), csh_aos.double_prime(tlo), std::numeric_limits<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), csh_aos.double_prime(thi), 400*std::numeric_limits<Real>::epsilon());
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>