From 5cbfdb554abbf5f7f7f1fe034b6345bc5335f39f Mon Sep 17 00:00:00 2001 From: Nick Date: Sun, 22 Mar 2020 11:22:36 -0400 Subject: [PATCH] Fix invalid read in cubic Hermite. [CI SKIP] --- doc/interpolators/cubic_hermite.qbk | 10 ++ .../math/interpolators/cubic_hermite.hpp | 17 +++- .../detail/cubic_hermite_detail.hpp | 29 +++++- test/cubic_hermite_test.cpp | 93 ++++++++++++++++--- test/daubechies_scaling_test.cpp | 8 ++ test/daubechies_wavelet_test.cpp | 21 +++++ 6 files changed, 164 insertions(+), 14 deletions(-) diff --git a/doc/interpolators/cubic_hermite.qbk b/doc/interpolators/cubic_hermite.qbk index 7b78e621e..347659c58 100644 --- a/doc/interpolators/cubic_hermite.qbk +++ b/doc/interpolators/cubic_hermite.qbk @@ -29,6 +29,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) void push_back(Real x, Real y, Real dydx); + std::pair domain() const; + friend std::ostream& operator<<(std::ostream & os, const cubic_hermite & m); }; @@ -42,6 +44,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) inline Real operator()(Real x) const; inline Real prime(Real x) const; + + std::pair domain() const; }; @@ -56,6 +60,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) inline Real operator()(Real x) const; inline Real prime(Real x) const; + + std::pair domain() const; }; } // namespaces @@ -124,6 +130,10 @@ For the "array of structs" version: auto ch = cardinal_cubic_hermite_aos(std::move(data), x0, dx); +For a quick sanity check, we can call `ch.domain()`: + + auto [x_min, x_max] = ch.domain(); + [heading Performance] Google benchmark was used to evaluate the performance. diff --git a/include/boost/math/interpolators/cubic_hermite.hpp b/include/boost/math/interpolators/cubic_hermite.hpp index f7a0bed07..c568510ec 100644 --- a/include/boost/math/interpolators/cubic_hermite.hpp +++ b/include/boost/math/interpolators/cubic_hermite.hpp @@ -39,11 +39,16 @@ public: impl_->push_back(x, y, dydx); } - int64_t bytes() + int64_t bytes() const { return impl_->bytes() + sizeof(impl_); } + std::pair domain() const + { + return impl_->domain(); + } + private: std::shared_ptr> impl_; }; @@ -78,6 +83,11 @@ public: return impl_->bytes() + sizeof(impl_); } + std::pair domain() const + { + return impl_->domain(); + } + private: std::shared_ptr> impl_; }; @@ -114,6 +124,11 @@ public: return impl_->bytes() + sizeof(impl_); } + std::pair domain() const + { + return impl_->domain(); + } + private: std::shared_ptr> impl_; }; diff --git a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp index 94c250ff3..6365e1ab6 100644 --- a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp @@ -155,6 +155,11 @@ public: return 3*x_.size()*sizeof(Real) + 3*sizeof(x_); } + std::pair domain() const + { + return {x_.front(), x_.back()}; + } + RandomAccessContainer x_; RandomAccessContainer y_; RandomAccessContainer dydx_; @@ -273,6 +278,12 @@ public: return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real); } + std::pair domain() const + { + Real xf = x0_ + (y_.size()-1)/inv_dx_; + return {x0_, xf}; + } + private: RandomAccessContainer y_; @@ -336,6 +347,12 @@ public: auto i = static_cast(ii); Real t = s - ii; + // If we had infinite precision, this would never happen. + // But we don't have infinite precision. + if (t == 0) + { + return dat_[i][0]; + } Real y0 = dat_[i][0]; Real y1 = dat_[i+1][0]; Real dy0 = dat_[i][1]; @@ -359,7 +376,7 @@ public: } if (x == xf) { - return dat_.back()[1]; + return dat_.back()[1]*inv_dx_; } return this->unchecked_prime(x); } @@ -371,6 +388,10 @@ public: Real ii = floor(s); auto i = static_cast(ii); Real t = s - ii; + if (t == 0) + { + return dat_[i][1]*inv_dx_; + } Real y0 = dat_[i][0]; Real dy0 = dat_[i][1]; Real y1 = dat_[i+1][0]; @@ -391,6 +412,12 @@ public: return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real); } + std::pair domain() const + { + Real xf = x0_ + (dat_.size()-1)/inv_dx_; + return {x0_, xf}; + } + private: RandomAccessContainer dat_; diff --git a/test/cubic_hermite_test.cpp b/test/cubic_hermite_test.cpp index 6d75c2883..4c9cfc003 100644 --- a/test/cubic_hermite_test.cpp +++ b/test/cubic_hermite_test.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #ifdef BOOST_HAS_FLOAT128 #include @@ -27,10 +28,11 @@ using boost::math::interpolators::cardinal_cubic_hermite_aos; template void test_constant() { - - std::vector x{0,1,2,3, 9, 22, 81}; + Real x0 = 0; + std::vector x{x0,1,2,3, 9, 22, 81}; std::vector y(x.size()); - for (auto & t : y) { + for (auto & t : y) + { t = 7; } @@ -40,9 +42,19 @@ void test_constant() auto dydx_copy = dydx; auto hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy)); - for (Real t = x[0]; t <= x.back(); t += 0.25) { - CHECK_ULP_CLOSE(Real(7), hermite_spline(t), 2); - CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(t), 2); + // Now check the boundaries: + Real tlo = x.front(); + Real thi = x.back(); + int samples = 5000; + int i = 0; + while (i++ < samples) + { + CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2); + CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2); + tlo = boost::math::nextafter(tlo, std::numeric_limits::max()); + thi = boost::math::nextafter(thi, std::numeric_limits::lowest()); } boost::circular_buffer x_buf(x.size()); @@ -226,10 +238,34 @@ void test_cardinal_constant() auto hermite_spline_aos = cardinal_cubic_hermite_aos(std::move(data), x0, dx); for (Real t = x0; t <= x0 + 24*dx; t += 0.25) { - CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(t), 2); - CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(t), 2); + if (!CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(t), 2)) { + std::cerr << " Wrong evaluation at t = " << t << "\n"; + } + if (!CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(t), 2)) { + std::cerr << " Wrong evaluation at t = " << t << "\n"; + } } - + + // Now check the boundaries: + Real tlo = x0; + Real thi = x0 + (25-1)*dx; + int samples = 5000; + int i = 0; + while (i++ < samples) + { + CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2); + CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2); + CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(tlo), 2); + CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(thi), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(tlo), 2); + CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(thi), 2); + + tlo = boost::math::nextafter(tlo, std::numeric_limits::max()); + thi = boost::math::nextafter(thi, std::numeric_limits::lowest()); + } + } @@ -278,6 +314,25 @@ void test_cardinal_linear() CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(t), 0); } + Real tlo = x0; + Real thi = x0 + (45-1)*dx; + int samples = 5000; + int i = 0; + while (i++ < samples) + { + CHECK_ULP_CLOSE(Real(tlo), hermite_spline(tlo), 2); + CHECK_ULP_CLOSE(Real(thi), hermite_spline(thi), 2); + CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(tlo), 2); + CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(thi), 2); + CHECK_ULP_CLOSE(Real(tlo), hermite_spline_aos(tlo), 2); + CHECK_ULP_CLOSE(Real(thi), hermite_spline_aos(thi), 2); + CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(tlo), 2); + CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(thi), 2); + + tlo = boost::math::nextafter(tlo, std::numeric_limits::max()); + thi = boost::math::nextafter(thi, std::numeric_limits::lowest()); + } + } @@ -318,6 +373,23 @@ void test_cardinal_quadratic() CHECK_ULP_CLOSE(t, saos.prime(t), 70); } + auto [tlo, thi] = s.domain(); + int samples = 5000; + int i = 0; + while (i++ < samples) + { + CHECK_ULP_CLOSE(Real(tlo*tlo/2), s(tlo), 3); + CHECK_ULP_CLOSE(Real(thi*thi/2), s(thi), 3); + CHECK_ULP_CLOSE(Real(tlo), s.prime(tlo), 3); + CHECK_ULP_CLOSE(Real(thi), s.prime(thi), 3); + CHECK_ULP_CLOSE(Real(tlo*tlo/2), saos(tlo), 3); + CHECK_ULP_CLOSE(Real(thi*thi/2), saos(thi), 3); + CHECK_ULP_CLOSE(Real(tlo), saos.prime(tlo), 3); + CHECK_ULP_CLOSE(Real(thi), saos.prime(thi), 3); + + tlo = boost::math::nextafter(tlo, std::numeric_limits::max()); + thi = boost::math::nextafter(thi, std::numeric_limits::lowest()); + } } @@ -359,8 +431,6 @@ int main() test_cardinal_quadratic(); test_cardinal_interpolation_condition(); - - test_constant(); test_linear(); test_quadratic(); @@ -370,7 +440,6 @@ int main() test_cardinal_quadratic(); test_cardinal_interpolation_condition(); - test_constant(); test_linear(); test_quadratic(); diff --git a/test/daubechies_scaling_test.cpp b/test/daubechies_scaling_test.cpp index 668cd609c..ceb43da88 100644 --- a/test/daubechies_scaling_test.cpp +++ b/test/daubechies_scaling_test.cpp @@ -427,6 +427,14 @@ void test_quadratures() { CHECK_ULP_CLOSE(Real(0), phi(xlo), 0); CHECK_ULP_CLOSE(Real(0), phi(xhi), 0); + if constexpr (p > 2) { + assert(abs(phi.prime(xlo)) <= 5); + assert(abs(phi.prime(xhi)) <= 5); + if constexpr (p > 5) { + assert(abs(phi.double_prime(xlo)) <= 5); + assert(abs(phi.double_prime(xhi)) <= 5); + } + } xlo = std::nextafter(xlo, std::numeric_limits::lowest()); xhi = std::nextafter(xhi, std::numeric_limits::max()); } diff --git a/test/daubechies_wavelet_test.cpp b/test/daubechies_wavelet_test.cpp index f1614d10e..384b02664 100644 --- a/test/daubechies_wavelet_test.cpp +++ b/test/daubechies_wavelet_test.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -81,6 +82,15 @@ void test_quadratures() { CHECK_ULP_CLOSE(Real(0), psi(xlo), 0); CHECK_ULP_CLOSE(Real(0), psi(xhi), 0); + if constexpr (p > 2) + { + CHECK_ULP_CLOSE(Real(0), psi.prime(xlo), 0); + CHECK_ULP_CLOSE(Real(0), psi.prime(xhi), 0); + if constexpr (p >= 6) { + CHECK_ULP_CLOSE(Real(0), psi.double_prime(xlo), 0); + CHECK_ULP_CLOSE(Real(0), psi.double_prime(xhi), 0); + } + } xlo = std::nextafter(xlo, std::numeric_limits::lowest()); xhi = std::nextafter(xhi, std::numeric_limits::max()); } @@ -88,8 +98,19 @@ void test_quadratures() xlo = a; xhi = b; for (int i = 0; i < samples; ++i) { + std::cout << std::setprecision(std::numeric_limits::max_digits10); assert(abs(psi(xlo)) <= 5); assert(abs(psi(xhi)) <= 5); + if constexpr (p > 2) + { + assert(abs(psi.prime(xlo)) <= 5); + assert(abs(psi.prime(xhi)) <= 5); + if constexpr (p >= 6) + { + assert(abs(psi.double_prime(xlo)) <= 5); + assert(abs(psi.double_prime(xhi)) <= 5); + } + } xlo = std::nextafter(xlo, std::numeric_limits::max()); xhi = std::nextafter(xhi, std::numeric_limits::lowest()); }