From fb26bc0cf19e3f645f4cac595f2ca1d6003d3588 Mon Sep 17 00:00:00 2001
From: Nick <nathompson7@protonmail.com>
Date: Sun, 22 Mar 2020 13:03:48 -0400
Subject: [PATCH] Fix septic hermite evaluation.

---
 .../detail/septic_hermite_detail.hpp          | 39 ++++++++++++++-
 .../math/interpolators/septic_hermite.hpp     | 15 ++++++
 test/septic_hermite_test.cpp                  | 48 +++++++++++++++++++
 3 files changed, 101 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
index 659cf019d..44843c6bb 100644
--- a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
+++ b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
@@ -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_;
diff --git a/include/boost/math/interpolators/septic_hermite.hpp b/include/boost/math/interpolators/septic_hermite.hpp
index 81d33291f..b5acbc764 100644
--- a/include/boost/math/interpolators/septic_hermite.hpp
+++ b/include/boost/math/interpolators/septic_hermite.hpp
@@ -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_;
 };
diff --git a/test/septic_hermite_test.cpp b/test/septic_hermite_test.cpp
index 6b7265327..944a3d0c4 100644
--- a/test/septic_hermite_test.cpp
+++ b/test/septic_hermite_test.cpp
@@ -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>