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 quintic hermite interpolation [CI SKIP]

Browse Source
This commit is contained in:
Nick
2020-03-22 12:07:05 -04:00
parent 5cbfdb554a
commit 2af33f2725
5 changed files with 119 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/interpolators/quintic_hermite.qbk
View File

@@ -25,6 +25,8 @@ public:
inline Real double_prime(Real x) const;
std::pair<Real, Real> domain() const;
friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m);
void push_back(Real x, Real y, Real dydx, Real d2ydx2);
@@ -41,6 +43,8 @@ public:
inline Real prime(Real x) const;
inline Real double_prime(Real x) const;
std::pair<Real, Real> domain() const;
};
template<class RandomAccessContainer>
@@ -56,6 +60,8 @@ public:
inline Real double_prime(Real x) const;
std::pair<Real, Real> domain() const;
}
``

45
include/boost/math/interpolators/detail/quintic_hermite_detail.hpp
View File

@@ -188,6 +188,11 @@ public:
return 4*x_.size()*sizeof(x_);
}
std::pair<Real, Real> domain() const
{
return {x_.front(), x_.back()};
}
private:
RandomAccessContainer x_;
RandomAccessContainer y_;
@@ -257,7 +262,10 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s - ii;
if (t == 0)
{
return y_[i];
}
Real y0 = y_[i];
Real y1 = y_[i+1];
Real dy0 = dy_[i];
@@ -301,7 +309,10 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s - ii;
if (t == 0)
{
return dy_[i]*inv_dx_;
}
Real y0 = y_[i];
Real y1 = y_[i+1];
Real dy0 = dy_[i];
@@ -341,7 +352,10 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s - ii;
if (t==0)
{
return d2y_[i]*2*inv_dx_*inv_dx_;
}
Real y0 = y_[i];
Real y1 = y_[i+1];
@@ -361,6 +375,12 @@ public:
return 3*y_.size()*sizeof(Real) + 2*sizeof(Real);
}
std::pair<Real, Real> domain() const
{
Real xf = x0_ + (y_.size()-1)/inv_dx_;
return {x0_, xf};
}
private:
RandomAccessContainer y_;
RandomAccessContainer dy_;
@@ -425,6 +445,10 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s - ii;
if (t == 0)
{
return data_[i][0];
}
Real y0 = data_[i][0];
Real dy0 = data_[i][1];
@@ -466,6 +490,10 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s - ii;
if (t == 0)
{
return data_[i][1]*inv_dx_;
}
Real y0 = data_[i][0];
@@ -507,8 +535,9 @@ public:
Real ii = floor(s);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s - ii;
if (t == 0) {
return data_[i][2]*2*inv_dx_*inv_dx_;
}
Real y0 = data_[i][0];
Real dy0 = data_[i][1];
Real d2y0 = data_[i][2];
@@ -527,6 +556,12 @@ public:
return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real);
}
std::pair<Real, Real> domain() const
{
Real xf = x0_ + (data_.size()-1)/inv_dx_;
return {x0_, xf};
}
private:
RandomAccessContainer data_;
Real x0_;

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

@@ -46,6 +46,11 @@ public:
return impl_->bytes() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::quintic_hermite_detail<RandomAccessContainer>> impl_;
};
@@ -76,6 +81,11 @@ public:
return impl_->bytes() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::cardinal_quintic_hermite_detail<RandomAccessContainer>> impl_;
};
@@ -108,6 +118,11 @@ public:
{
return impl_->bytes() + sizeof(impl_);
}
std::pair<Real, Real> domain() const
{
return impl_->domain();
}
private:
std::shared_ptr<detail::cardinal_quintic_hermite_detail_aos<RandomAccessContainer>> impl_;
};

2
test/cubic_hermite_test.cpp
View File

@@ -246,7 +246,7 @@ void test_cardinal_constant()
}
}
// Now check the boundaries:
// Now check the boundaries:
Real tlo = x0;
Real thi = x0 + (25-1)*dx;
int samples = 5000;

57
test/quintic_hermite_test.cpp
View File

@@ -13,6 +13,7 @@
#include <boost/random/uniform_real.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/math/interpolators/quintic_hermite.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/circular_buffer.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
@@ -287,6 +288,25 @@ void test_cardinal_constant()
CHECK_ULP_CLOSE(Real(0), qh_aos.prime(t), 24);
CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 24);
}
// Now check the boundaries:
auto [tlo, thi] = qh.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(7), qh(tlo), 2);
CHECK_ULP_CLOSE(Real(7), qh(thi), 2);
CHECK_ULP_CLOSE(Real(7), qh_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(7), qh_aos(thi), 2);
CHECK_ULP_CLOSE(Real(0), qh.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), qh.prime(thi), 2);
CHECK_ULP_CLOSE(Real(0), qh_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(0), qh_aos.prime(thi), 2);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
@@ -322,6 +342,24 @@ void test_cardinal_linear()
CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 2);
}
// Now check the boundaries:
auto [tlo, thi] = qh.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(Real(tlo), qh(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), qh(thi), 2);
CHECK_ULP_CLOSE(Real(tlo), qh_aos(tlo), 2);
CHECK_ULP_CLOSE(Real(thi), qh_aos(thi), 2);
CHECK_ULP_CLOSE(Real(1), qh.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), qh.prime(thi), 128);
CHECK_ULP_CLOSE(Real(1), qh_aos.prime(tlo), 2);
CHECK_ULP_CLOSE(Real(1), qh_aos.prime(thi), 128);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>
@@ -366,6 +404,25 @@ void test_cardinal_quadratic()
CHECK_ULP_CLOSE(t, qh_aos.prime(t), 12);
CHECK_ULP_CLOSE(Real(1), qh_aos.double_prime(t), 64);
}
// Now check the boundaries:
auto [tlo, thi] = qh.domain();
int samples = 5000;
int i = 0;
while (i++ < samples)
{
CHECK_ULP_CLOSE(tlo*tlo/2, qh(tlo), 16);
CHECK_ULP_CLOSE(thi*thi/2, qh(thi), 16);
CHECK_ULP_CLOSE(tlo*tlo/2, qh_aos(tlo), 16);
CHECK_ULP_CLOSE(thi*thi/2, qh_aos(thi), 16);
CHECK_ULP_CLOSE(tlo, qh.prime(tlo), 16);
CHECK_ULP_CLOSE(thi, qh.prime(thi), 64);
CHECK_ULP_CLOSE(tlo, qh_aos.prime(tlo), 16);
CHECK_ULP_CLOSE(thi, qh_aos.prime(thi), 64);
tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
}
}
template<typename Real>