From de67bcb45c7fa4b17ee5fb2b5e9ca42cddf8c80e Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 16 Jan 2020 16:31:32 -0500 Subject: [PATCH 1/5] Quintic Hermite interpolation [CI SKIP] --- doc/interpolators/quintic_hermite.qbk | 81 ++++++++++ doc/math.qbk | 1 + .../detail/quintic_hermite_detail.hpp | 131 ++++++++++++++++ .../math/interpolators/quintic_hermite.hpp | 41 +++++ test/quintic_hermite_test.cpp | 144 ++++++++++++++++++ 5 files changed, 398 insertions(+) create mode 100644 doc/interpolators/quintic_hermite.qbk create mode 100644 include/boost/math/interpolators/detail/quintic_hermite_detail.hpp create mode 100644 include/boost/math/interpolators/quintic_hermite.hpp create mode 100644 test/quintic_hermite_test.cpp diff --git a/doc/interpolators/quintic_hermite.qbk b/doc/interpolators/quintic_hermite.qbk new file mode 100644 index 000000000..7973f2c42 --- /dev/null +++ b/doc/interpolators/quintic_hermite.qbk @@ -0,0 +1,81 @@ +[/ +Copyright (c) 2020 Nick Thompson +Use, modification and distribution are subject to the +Boost Software License, Version 1.0. (See accompanying file +LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +] + +[section:quintic_hermite Quintic Hermite interpolation] + +[heading Synopsis] +`` +#include + +namespace boost::math::interpolators { + +template +class quintic_hermite { +public: + using Real = typename RandomAccessContainer::value_type; + quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) + + Real operator()(Real x) const; + + Real prime(Real x) const; + + friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m); + + void push_back(Real x, Real y, Real dydx, Real d2ydx2); +}; +} +`` + + +[heading Quintic Hermite Interpolation] + +The quintic Hermite interpolator takes a list of possibly non-uniformly spaced abscissas, ordinates, and their velocities and accelerations which are used to construct a quintic interpolating polynomial between segments. +The interpolant is /C/[super 2] and its evaluation has [bigo](log(/N/)) complexity. +An example usage is as follows: + + std::vector x{1,2,3, 4, 5, 6}; + std::vector y{7,8,9,10,11,12}; + std::vector dydx{1,1,1,1,1,1}; + std::vector d2ydx2{0,0,0,0,0,0}; + + using boost::math::interpolators::quintic_hermite; + auto spline = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); + // evaluate at a point: + double z = spline(3.4); + // evaluate derivative at a point: + double zprime = spline.prime(3.4); + +Periodically, it is helpful to see what data the interpolator has. +This can be achieved via + + std::cout << spline << "\n"; + +Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads. +(The call operator and `.prime()` are threadsafe.) + +One unique aspect of this interpolator is that it can be updated in constant time. +Hence we can use `boost::circular_buffer` to do real-time interpolation. + + + #include + ... + boost::circular_buffer initial_x{1,2,3,4}; + boost::circular_buffer initial_y{4,5,6,7}; + boost::circular_buffer initial_dydx{1,1,1,1}; + boost::circular_buffer initial_d2ydx2{0,0,0,0}; + auto circular_akima = quintic_hermite(std::move(initial_x), std::move(initial_y), std::move(initial_dydx), std::move(initial_d2ydx2)); + // interpolate via call operation: + double y = circular_akima(3.5); + // add new data: + circular_akima.push_back(5, 8, 1, 0); + // interpolat at 4.5: + y = circular_akima(4.5); + + + +[endsect] +[/section:makima] diff --git a/doc/math.qbk b/doc/math.qbk index b5c29b8c8..91bcdcdb7 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -714,6 +714,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include interpolators/cardinal_trigonometric.qbk] [include interpolators/makima.qbk] [include interpolators/pchip.qbk] +[include interpolators/quintic_hermite.qbk] [endmathpart] [mathpart quadrature Quadrature and Differentiation] diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp new file mode 100644 index 000000000..2680cca84 --- /dev/null +++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp @@ -0,0 +1,131 @@ +#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP +#include +#include + +namespace boost::math::interpolators::detail { + +template +class quintic_hermite_detail { +public: + using Real = typename RandomAccessContainer::value_type; + quintic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)} + { + if (x_.size() != y_.size()) { + throw std::domain_error("Number of abscissas must = number of ordinates."); + } + if (x_.size() != dydx_.size()) { + throw std::domain_error("Numbers of derivatives must = number of abscissas."); + } + if (x_.size() != d2ydx2_.size()) { + throw std::domain_error("Number of second derivatives must equal number of abscissas."); + } + if (x_.size() < 2) { + throw std::domain_error("At least 2 abscissas are required."); + } + Real x0 = x_[0]; + for (decltype(x_.size()) i = 1; i < x_.size(); ++i) { + Real x1 = x_[i]; + if (x1 <= x0) + { + throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}"); + } + x0 = x1; + } + } + + void push_back(Real x, Real y, Real dydx, Real d2ydx2) { + using std::abs; + using std::isnan; + if (x <= x_.back()) { + throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); + } + x_.push_back(x); + y_.push_back(y); + dydx_.push_back(dydx); + d2ydx2_.push_back(d2ydx2); + } + + Real operator()(Real x) const { + if (x < x_[0] || x > x_.back()) { + std::ostringstream oss; + oss.precision(std::numeric_limits::digits10+3); + oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" + << x_[0] << ", " << x_.back() << "]"; + throw std::domain_error(oss.str()); + } + // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work. + // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf. + if (x == x_.back()) { + return y_.back(); + } + + auto it = std::upper_bound(x_.begin(), x_.end(), x); + auto i = std::distance(x_.begin(), it) -1; + Real x0 = *(it-1); + Real x1 = *it; + Real y0 = y_[i]; + Real y1 = y_[i+1]; + Real v0 = dydx_[i]; + Real v1 = dydx_[i+1]; + Real a0 = d2ydx2_[i]; + Real a1 = d2ydx2_[i+1]; + + Real dx = (x1-x0); + Real t = (x-x0)/dx; + + // See the 'Basis functions' section of: + // https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf + // Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs + Real y = (1- t*t*t*(10 - 15*t + 6*t*t))*y0 + (t-6*t*t*t + 8*t*t*t*t -3*t*t*t*t*t)*v0 + (t*t-3*t*t*t +3*t*t*t*t-t*t*t*t*t)*a0/2; + y += t*t*t*(1 - 2*t + t*t)*a1/2 + (-4*t*t*t + 7*t*t*t*t -3*t*t*t*t*t)*v1 + (10*t*t*t - 15*t*t*t*t + 6*t*t*t*t*t)*y1; + return y; + } + + Real prime(Real x) const { + if (x < x_[0] || x > x_.back()) { + std::ostringstream oss; + oss.precision(std::numeric_limits::digits10+3); + oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" + << x_[0] << ", " << x_.back() << "]"; + throw std::domain_error(oss.str()); + } + if (x == x_.back()) { + return dydx_.back(); + } + + auto it = std::upper_bound(x_.begin(), x_.end(), x); + auto i = std::distance(x_.begin(), it) -1; + Real x0 = *(it-1); + Real x1 = *it; + Real s0 = dydx_[i]; + Real s1 = dydx_[i+1]; + + // Ridiculous linear interpolation. Fine for now: + Real numerator = s0*(x1-x) + s1*(x-x0); + Real denominator = x1 - x0; + return numerator/denominator; + } + + + friend std::ostream& operator<<(std::ostream & os, const quintic_hermite_detail & m) + { + os << "(x,y,y') = {"; + for (size_t i = 0; i < m.x_.size() - 1; ++i) { + os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] << "), "; + } + auto n = m.x_.size()-1; + os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << ")}"; + return os; + } + + + +private: + RandomAccessContainer x_; + RandomAccessContainer y_; + RandomAccessContainer dydx_; + RandomAccessContainer d2ydx2_; +}; +} +#endif \ No newline at end of file diff --git a/include/boost/math/interpolators/quintic_hermite.hpp b/include/boost/math/interpolators/quintic_hermite.hpp new file mode 100644 index 000000000..0111df9a1 --- /dev/null +++ b/include/boost/math/interpolators/quintic_hermite.hpp @@ -0,0 +1,41 @@ +#ifndef BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP +#define BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP +#include +#include +#include +#include + +namespace boost::math::interpolators { + +template +class quintic_hermite { +public: + using Real = typename RandomAccessContainer::value_type; + quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) + : impl_(std::make_shared>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2))) + {} + + Real operator()(Real x) const { + return impl_->operator()(x); + } + + Real prime(Real x) const { + return impl_->prime(x); + } + + friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m) + { + os << *m.impl_; + return os; + } + + void push_back(Real x, Real y, Real dydx, Real d2ydx2) { + impl_->push_back(x, y, dydx, d2ydx2); + } + + +private: + std::shared_ptr> impl_; +}; +} +#endif \ No newline at end of file diff --git a/test/quintic_hermite_test.cpp b/test/quintic_hermite_test.cpp new file mode 100644 index 000000000..4e5e9dd26 --- /dev/null +++ b/test/quintic_hermite_test.cpp @@ -0,0 +1,144 @@ +/* + * Copyright Nick Thompson, 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#include "math_unit_test.hpp" +#include +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + + +using boost::math::interpolators::quintic_hermite; + +template +void test_constant() +{ + + std::vector x{0,1,2,3, 9, 22, 81}; + std::vector y(x.size()); + std::vector dydx(x.size(), 0); + std::vector d2ydx2(x.size(), 0); + for (auto & t : y) { + t = 7; + } + + auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); + + for (Real t = 0; t <= 81; t += 0.25) { + CHECK_ULP_CLOSE(Real(7), qh(t), 24); + CHECK_ULP_CLOSE(Real(0), qh.prime(t), 24); + } +} + + +template +void test_linear() +{ + + 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); + + auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); + + for (Real t = 0; t <= 9; t += 0.25) { + CHECK_ULP_CLOSE(Real(t), qh(t), 2); + CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2); + } +} + +template +void test_quadratic() +{ + + std::vector x{0,1,2,3, 4,5,6,7,8,9}; + std::vector y(x.size()); + for (size_t i = 0; i < y.size(); ++i) + { + y[i] = x[i]*x[i]/2; + } + + std::vector dydx(x.size()); + for (size_t i = 0; i < y.size(); ++i) { + dydx[i] = x[i]; + } + + std::vector d2ydx2(x.size(), 1); + + auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); + + for (Real t = 0; t <= 9; t += 0.125) { + CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 2); + CHECK_ULP_CLOSE(t, qh.prime(t), 2); + } +} + +template +void test_interpolation_condition() +{ + for (size_t n = 4; n < 50; ++n) { + std::vector x(n); + std::vector y(n); + std::vector dydx(n); + std::vector d2ydx2(n); + std::default_random_engine rd; + std::uniform_real_distribution dis(0,1); + Real x0 = dis(rd); + x[0] = x0; + y[0] = dis(rd); + for (size_t i = 1; i < n; ++i) { + x[i] = x[i-1] + dis(rd); + y[i] = dis(rd); + dydx[i] = dis(rd); + d2ydx2[i] = dis(rd); + } + + auto x_copy = x; + auto y_copy = y; + auto dydx_copy = dydx; + auto d2ydx2_copy = d2ydx2; + auto s = quintic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy)); + //std::cout << "s = " << s << "\n"; + for (size_t i = 0; i < x.size(); ++i) { + CHECK_ULP_CLOSE(y[i], s(x[i]), 2); + CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2); + } + } +} + + +int main() +{ + test_constant(); + test_linear(); + test_quadratic(); + test_interpolation_condition(); + + test_constant(); + test_linear(); + test_quadratic(); + test_interpolation_condition(); + + test_constant(); + test_linear(); + test_quadratic(); + test_interpolation_condition(); + +#ifdef BOOST_HAS_FLOAT128 + test_constant(); + test_linear(); + test_quadratic(); +#endif + + return boost::math::test::report_errors(); +} From 0242d1ac6b1048f422cd245777a68f82b33cdbe7 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 16 Jan 2020 18:08:34 -0500 Subject: [PATCH 2/5] Quintic Hermite: Fix #include, fix docs, add test to Jamfile. --- doc/interpolators/quintic_hermite.qbk | 46 +++++++++++++++++-- .../detail/quintic_hermite_detail.hpp | 1 + test/Jamfile.v2 | 1 + 3 files changed, 45 insertions(+), 3 deletions(-) diff --git a/doc/interpolators/quintic_hermite.qbk b/doc/interpolators/quintic_hermite.qbk index 7973f2c42..17ff0d017 100644 --- a/doc/interpolators/quintic_hermite.qbk +++ b/doc/interpolators/quintic_hermite.qbk @@ -57,7 +57,7 @@ This can be achieved via Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads. (The call operator and `.prime()` are threadsafe.) -One unique aspect of this interpolator is that it can be updated in constant time. +The interpolator can be updated in constant time. Hence we can use `boost::circular_buffer` to do real-time interpolation. @@ -72,10 +72,50 @@ Hence we can use `boost::circular_buffer` to do real-time interpolation. double y = circular_akima(3.5); // add new data: circular_akima.push_back(5, 8, 1, 0); - // interpolat at 4.5: + // interpolate at 4.5: y = circular_akima(4.5); +[heading Complexity and Performance] +The following google benchmark demonstrates the cost of the call operator for this interpolator: + +``` +Running ./bench.x +Run on (16 X 4300 MHz CPU s) +CPU Caches: + L1 Data 32K (x8) + L1 Instruction 32K (x8) + L2 Unified 1024K (x8) + L3 Unified 11264K (x1) +Load Average: 0.92, 0.64, 0.35 +***WARNING*** CPU scaling is enabled, the benchmark real time measurements may be noisy and will incur extra overhead. +----------------------------------------------- +Benchmark Time +----------------------------------------------- +BMQuinticHermite/8 10.5 ns +BMQuinticHermite/16 11.0 ns +BMQuinticHermite/32 11.5 ns +BMQuinticHermite/64 12.1 ns +BMQuinticHermite/128 12.8 ns +BMQuinticHermite/256 13.3 ns +BMQuinticHermite/512 13.9 ns +BMQuinticHermite/1024 14.6 ns +BMQuinticHermite/2048 15.4 ns +BMQuinticHermite/4096 16.3 ns +BMQuinticHermite/8192 17.5 ns +BMQuinticHermite/16384 18.8 ns +BMQuinticHermite/32768 20.0 ns +BMQuinticHermite/65536 21.1 ns +BMQuinticHermite/131072 22.1 ns +BMQuinticHermite/262144 23.2 ns +BMQuinticHermite/524288 24.2 ns +BMQuinticHermite/1048576 25.9 ns +BMQuinticHermite/2097152 27.5 ns +BMQuinticHermite/4194304 29.0 ns +BMQuinticHermite/8388608 30.2 ns +BMQuinticHermite/16777216 31.7 ns +BMQuinticHermite_BigO 1.35 lgN +``` [endsect] -[/section:makima] +[/section:quintic_hermite] diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp index 2680cca84..ed9e97bc9 100644 --- a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp @@ -2,6 +2,7 @@ #define BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP #include #include +#include namespace boost::math::interpolators::detail { diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index a1ce0e7da..92129c6d9 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -947,6 +947,7 @@ test-suite misc : [ run cardinal_quintic_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run makima_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run pchip_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run quintic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=2 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_2 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ] From 32f7e35012a93b1f70ef506334e355e9a7613042 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 16 Jan 2020 18:09:58 -0500 Subject: [PATCH 3/5] Quintic Hermite: Remove irrelevant information from docs [CI SKIP] --- doc/interpolators/quintic_hermite.qbk | 2 -- 1 file changed, 2 deletions(-) diff --git a/doc/interpolators/quintic_hermite.qbk b/doc/interpolators/quintic_hermite.qbk index 17ff0d017..958a8a92e 100644 --- a/doc/interpolators/quintic_hermite.qbk +++ b/doc/interpolators/quintic_hermite.qbk @@ -80,7 +80,6 @@ Hence we can use `boost::circular_buffer` to do real-time interpolation. The following google benchmark demonstrates the cost of the call operator for this interpolator: ``` -Running ./bench.x Run on (16 X 4300 MHz CPU s) CPU Caches: L1 Data 32K (x8) @@ -88,7 +87,6 @@ CPU Caches: L2 Unified 1024K (x8) L3 Unified 11264K (x1) Load Average: 0.92, 0.64, 0.35 -***WARNING*** CPU scaling is enabled, the benchmark real time measurements may be noisy and will incur extra overhead. ----------------------------------------------- Benchmark Time ----------------------------------------------- From 408392dff03cacda69c7d34c594b05233b088ba0 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Fri, 17 Jan 2020 17:39:48 -0500 Subject: [PATCH 4/5] Update timing [CI SKIP] --- doc/interpolators/quintic_hermite.qbk | 46 +++++++++---------- .../detail/quintic_hermite_detail.hpp | 4 +- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/doc/interpolators/quintic_hermite.qbk b/doc/interpolators/quintic_hermite.qbk index 958a8a92e..1f008fbed 100644 --- a/doc/interpolators/quintic_hermite.qbk +++ b/doc/interpolators/quintic_hermite.qbk @@ -90,29 +90,29 @@ Load Average: 0.92, 0.64, 0.35 ----------------------------------------------- Benchmark Time ----------------------------------------------- -BMQuinticHermite/8 10.5 ns -BMQuinticHermite/16 11.0 ns -BMQuinticHermite/32 11.5 ns -BMQuinticHermite/64 12.1 ns -BMQuinticHermite/128 12.8 ns -BMQuinticHermite/256 13.3 ns -BMQuinticHermite/512 13.9 ns -BMQuinticHermite/1024 14.6 ns -BMQuinticHermite/2048 15.4 ns -BMQuinticHermite/4096 16.3 ns -BMQuinticHermite/8192 17.5 ns -BMQuinticHermite/16384 18.8 ns -BMQuinticHermite/32768 20.0 ns -BMQuinticHermite/65536 21.1 ns -BMQuinticHermite/131072 22.1 ns -BMQuinticHermite/262144 23.2 ns -BMQuinticHermite/524288 24.2 ns -BMQuinticHermite/1048576 25.9 ns -BMQuinticHermite/2097152 27.5 ns -BMQuinticHermite/4194304 29.0 ns -BMQuinticHermite/8388608 30.2 ns -BMQuinticHermite/16777216 31.7 ns -BMQuinticHermite_BigO 1.35 lgN +BMQuinticHermite/8 8.14 ns +BMQuinticHermite/16 8.69 ns +BMQuinticHermite/32 9.42 ns +BMQuinticHermite/64 9.90 ns +BMQuinticHermite/128 10.4 ns +BMQuinticHermite/256 10.9 ns +BMQuinticHermite/512 11.6 ns +BMQuinticHermite/1024 12.3 ns +BMQuinticHermite/2048 12.8 ns +BMQuinticHermite/4096 13.6 ns +BMQuinticHermite/8192 14.6 ns +BMQuinticHermite/16384 15.5 ns +BMQuinticHermite/32768 17.4 ns +BMQuinticHermite/65536 18.5 ns +BMQuinticHermite/131072 18.8 ns +BMQuinticHermite/262144 19.8 ns +BMQuinticHermite/524288 20.5 ns +BMQuinticHermite/1048576 21.6 ns +BMQuinticHermite/2097152 22.5 ns +BMQuinticHermite/4194304 24.2 ns +BMQuinticHermite/8388608 26.6 ns +BMQuinticHermite/16777216 26.7 ns +BMQuinticHermite_BigO 1.14 lgN ``` [endsect] diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp index ed9e97bc9..9d3cfddc6 100644 --- a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp @@ -78,8 +78,8 @@ public: // See the 'Basis functions' section of: // https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf // Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs - Real y = (1- t*t*t*(10 - 15*t + 6*t*t))*y0 + (t-6*t*t*t + 8*t*t*t*t -3*t*t*t*t*t)*v0 + (t*t-3*t*t*t +3*t*t*t*t-t*t*t*t*t)*a0/2; - y += t*t*t*(1 - 2*t + t*t)*a1/2 + (-4*t*t*t + 7*t*t*t*t -3*t*t*t*t*t)*v1 + (10*t*t*t - 15*t*t*t*t + 6*t*t*t*t*t)*y1; + Real y = (1- t*t*t*(10 - 15*t + 6*t*t))*y0 + t*(1-6*t*t + 8*t*t*t -3*t*t*t*t)*v0 + t*t*(1-3*t +3*t*t-t*t*t)*a0/2; + y += t*t*t*((1 - 2*t + t*t)*a1/2 + (-4 + 7*t -3*t*t)*v1 + (10 - 15*t + 6*t*t)*y1); return y; } From f1a90ae686b8ee14b9aa57d15a2664ffbe16cb22 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Tue, 21 Jan 2020 10:01:04 -0500 Subject: [PATCH 5/5] Quintic Hermite interpolation: Add test of cubic polynomials, fully Hornerize the basis functions. --- .../detail/quintic_hermite_detail.hpp | 7 ++-- test/quintic_hermite_test.cpp | 33 +++++++++++++++++++ 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp index 9d3cfddc6..bb313f820 100644 --- a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp @@ -3,6 +3,7 @@ #include #include #include +#include namespace boost::math::interpolators::detail { @@ -78,8 +79,10 @@ public: // See the 'Basis functions' section of: // https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf // Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs - Real y = (1- t*t*t*(10 - 15*t + 6*t*t))*y0 + t*(1-6*t*t + 8*t*t*t -3*t*t*t*t)*v0 + t*t*(1-3*t +3*t*t-t*t*t)*a0/2; - y += t*t*t*((1 - 2*t + t*t)*a1/2 + (-4 + 7*t -3*t*t)*v1 + (10 - 15*t + 6*t*t)*y1); + Real y = (1- t*t*t*(10 + t*(-15 + 6*t)))*y0; + y += t*(1+ t*t*(-6 + t*(8 -3*t)))*v0; + y += t*t*(1 + t*(-3 + t*(3-t)))*a0/2; + y += t*t*t*((1 + t*(-2 + t))*a1/2 + (-4 + t*(7 -3*t))*v1 + (10 + t*(-15 + 6*t))*y1); return y; } diff --git a/test/quintic_hermite_test.cpp b/test/quintic_hermite_test.cpp index 4e5e9dd26..befafbb44 100644 --- a/test/quintic_hermite_test.cpp +++ b/test/quintic_hermite_test.cpp @@ -83,6 +83,35 @@ void test_quadratic() } } +template +void test_cubic() +{ + + std::vector x{0,1,2,3, 4,5,6,7,8,9}; + std::vector y(x.size()); + for (size_t i = 0; i < y.size(); ++i) + { + y[i] = x[i]*x[i]*x[i]/6; + } + + std::vector dydx(x.size()); + for (size_t i = 0; i < y.size(); ++i) { + dydx[i] = x[i]*x[i]/2; + } + + std::vector d2ydx2(x.size()); + for (size_t i = 0; i < y.size(); ++i) { + d2ydx2[i] = x[i]; + } + + auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); + + for (Real t = 0; t <= 9; t += 0.125) { + CHECK_ULP_CLOSE(Real(t*t*t)/6, qh(t), 5); + //CHECK_ULP_CLOSE(t, qh.prime(t), 2); + } +} + template void test_interpolation_condition() { @@ -122,22 +151,26 @@ int main() test_constant(); test_linear(); test_quadratic(); + test_cubic(); test_interpolation_condition(); test_constant(); test_linear(); test_quadratic(); + test_cubic(); test_interpolation_condition(); test_constant(); test_linear(); test_quadratic(); + test_cubic(); test_interpolation_condition(); #ifdef BOOST_HAS_FLOAT128 test_constant(); test_linear(); test_quadratic(); + test_cubic(); #endif return boost::math::test::report_errors();