diff --git a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp index 8a8787980..122c2ae38 100644 --- a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp @@ -233,6 +233,28 @@ public: Real ii = floor(s3); auto i = static_cast(ii); Real t = s3 - ii; + // See: + // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m + Real t2 = t*t; + Real t3 = t2*t; + Real t4 = t3*t; + + Real s = t4*(-35 + t*(84 + t*(-70 + 20*t))); + Real z4 = -s; + Real z0 = s + 1; + //Real z1 = t*(10*t6 - 36*t5 + 45*t4 - 20*t3 + 1); + Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t)))); + //Real z2 = t2*(4*t4*t - 15*t4 + 20*t3 - 10*t2 + 1); + Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t)))); + //Real z3 = t3*(t4 - 4*t3 + 6*t2 - 4*t + 1); + Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t)))); + + //Real z5 = t4*(10*t3 - 34*t2 + 39*t - 15); + Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t))); + //Real z6 = t4*(-4*t3 + 13*t2 - 14*t + 5); + Real z6 = t4*(5 + t*(-14 + t*(13-4*t))); + //Real z7 = t4*(t3 - 3*t2 + 3*t - 1); + Real z7 = t4*(-1 + t*(3+t*(-3+t))); Real y0 = y_[i]; Real y1 = y_[i+1]; @@ -246,29 +268,6 @@ public: Real j0 = d3y_[i]; Real j1 = d3y_[i+1]; - // See: - // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m - Real t2 = t*t; - Real t3 = t2*t; - Real t4 = t3*t; - - Real s = t4*(-35 + t*(84 + t*(-70 + 20*t))); - Real z4 = -s; - Real z0 = s + 1; - //Real z1 = t*(10*t6 - 36*t5 + 45*t4 - 20*t3 + 1); - Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t)))); - //Real z2 = t2*(4*t5 - 15*t4 + 20*t3 - 10*t2 + 1); - Real z2 = t2*(1+ t2*(4*t3 - 15*t + 20*t - 10)); - //Real z3 = t3*(t4 - 4*t3 + 6*t2 - 4*t + 1); - Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t)))); - - //Real z5 = t4*(10*t3 - 34*t2 + 39*t - 15); - Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t))); - //Real z6 = t4*(-4*t3 + 13*t2 - 14*t + 5); - Real z6 = t4*(5 + t*(-14 + t*(13-4*t))); - //Real z7 = t4*(t3 - 3*t2 + 3*t - 1); - Real z7 = t4*(-1 + t*(3+t*(-3+t))); - return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1; } @@ -363,20 +362,6 @@ public: Real ii = floor(s3); auto i = static_cast(ii); Real t = s3 - ii; - - - Real y0 = data_[i][0]; - Real y1 = data_[i+1][0]; - // Velocity: - Real dy0 = data_[i][1]; - Real dy1 = data_[i+1][1]; - // Acceleration: - Real a0 = data_[i][2]; - Real a1 = data_[i+1][2]; - // Jerk: - Real j0 = data_[i][3]; - Real j1 = data_[i+1][3]; - Real t2 = t*t; Real t3 = t2*t; Real t4 = t3*t; @@ -387,7 +372,7 @@ public: //Real z1 = t*(10*t6 - 36*t5 + 45*t4 - 20*t3 + 1); Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t)))); //Real z2 = t2*(4*t5 - 15*t4 + 20*t3 - 10*t2 + 1); - Real z2 = t2*(1+ t2*(4*t3 - 15*t + 20*t - 10)); + Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t)))); //Real z3 = t3*(t4 - 4*t3 + 6*t2 - 4*t + 1); Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t)))); @@ -398,6 +383,15 @@ public: //Real z7 = t4*(t3 - 3*t2 + 3*t - 1); Real z7 = t4*(-1 + t*(3+t*(-3+t))); + Real y0 = data_[i][0]; + Real dy0 = data_[i][1]; + Real a0 = data_[i][2]; + Real j0 = data_[i][3]; + Real y1 = data_[i+1][0]; + Real dy1 = data_[i+1][1]; + Real a1 = data_[i+1][2]; + Real j1 = data_[i+1][3]; + return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1; } diff --git a/include/boost/math/interpolators/septic_hermite.hpp b/include/boost/math/interpolators/septic_hermite.hpp index 54da49a02..77ebe3024 100644 --- a/include/boost/math/interpolators/septic_hermite.hpp +++ b/include/boost/math/interpolators/septic_hermite.hpp @@ -57,5 +57,27 @@ private: std::shared_ptr> impl_; }; + +template +class cardinal_septic_hermite_aos { +public: + using Point = typename RandomAccessContainer::value_type; + using Real = typename Point::value_type; + cardinal_septic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx) + : impl_(std::make_shared>(std::move(data), x0, dx)) + {} + + Real operator()(Real x) const { + return impl_->operator()(x); + } + + Real prime(Real x) const { + return impl_->prime(x); + } + +private: + std::shared_ptr> impl_; +}; + } #endif \ No newline at end of file diff --git a/test/septic_hermite_test.cpp b/test/septic_hermite_test.cpp index ed9772328..6a7c4c7ea 100644 --- a/test/septic_hermite_test.cpp +++ b/test/septic_hermite_test.cpp @@ -19,6 +19,7 @@ using boost::multiprecision::float128; using boost::math::interpolators::septic_hermite; using boost::math::interpolators::cardinal_septic_hermite; +using boost::math::interpolators::cardinal_septic_hermite_aos; template void test_constant() @@ -35,7 +36,8 @@ void test_constant() auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3)); - for (Real t = 0; t <= 81; t += 0.25) { + for (Real t = 0; t <= 81; t += 0.25) + { CHECK_ULP_CLOSE(Real(7), sh(t), 24); CHECK_ULP_CLOSE(Real(0), sh.prime(t), 24); } @@ -52,13 +54,27 @@ void test_constant() //CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24); } + std::vector> data(128); + for (size_t i = 0; i < data.size(); ++i) + { + data[i][0] = 7; + data[i][1] = 0; + data[i][2] = 0; + data[i][3] = 0; + } + auto csh_aos = cardinal_septic_hermite_aos(std::move(data), x0, dx); + for (Real t = x0; t <= 127; t += 0.25) { + CHECK_ULP_CLOSE(Real(7), csh_aos(t), 24); + //CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24); + } + } template void test_linear() { - std::vector x{0,1,2,3, 4,5,6,7,8,9}; + std::vector x{0,1,2,3,4,5,6,7,8,9}; std::vector y = x; std::vector dydx(x.size(), 1); std::vector d2ydx2(x.size(), 0); @@ -88,9 +104,39 @@ void test_linear() sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3)); - for (Real t = xmin; t <= xmax; t += 0.125) { + for (Real t = xmin; t <= xmax; t += 0.125) + { CHECK_ULP_CLOSE(t, sh(t), 15); } + + Real x0 = 0; + Real dx = 1; + y.resize(10); + dydx.resize(10, 1); + d2ydx2.resize(10, 0); + d3ydx3.resize(10, 0); + for (size_t i = 0; i < y.size(); ++i) { + y[i] = i; + } + auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx); + for (Real t = 0; t <= 9; t += 0.125) + { + CHECK_ULP_CLOSE(t, csh(t), 15); + } + + std::vector> data(10); + for (size_t i = 0; i < data.size(); ++i) + { + data[i][0] = i; + data[i][1] = 1; + data[i][2] = 0; + data[i][3] = 0; + } + auto csh_aos = cardinal_septic_hermite_aos(std::move(data), x0, dx); + for (Real t = 0; t <= 9; t += 0.125) + { + CHECK_ULP_CLOSE(t, csh_aos(t), 15); + } } template @@ -148,8 +194,48 @@ void test_quadratic() CHECK_ULP_CLOSE(t*t/2, sh(t), 24); //CHECK_ULP_CLOSE(t, qh.prime(t), 36); } + + y.resize(10); + for (size_t i = 0; i < y.size(); ++i) + { + y[i] = i*i/Real(2); + } + + dydx.resize(y.size()); + for (size_t i = 0; i < y.size(); ++i) + { + dydx[i] = i; + } + + d2ydx2.resize(y.size(), 1); + d3ydx3.resize(y.size(), 0); + + Real x0 = 0; + Real dx = 1; + auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx); + for (Real t = x0; t <= 9; t += 0.125) + { + CHECK_ULP_CLOSE(t*t/2, csh(t), 24); + } + + std::vector> data(10); + for (size_t i = 0; i < data.size(); ++i) + { + data[i][0] = i*i/Real(2); + data[i][1] = i; + data[i][2] = 1; + data[i][3] = 0; + } + auto csh_aos = cardinal_septic_hermite_aos(std::move(data), x0, dx); + for (Real t = x0; t <= 9; t += 0.125) + { + CHECK_ULP_CLOSE(t*t/2, csh_aos(t), 24); + } + } + + template void test_cubic() {