Fix bug in septic_hermite; add unit tests [CI SKIP]

include/boost/math/interpolators/detail/septic_hermite_detail.hpp

@@ -233,6 +233,28 @@ public:
         Real ii = floor(s3);
         auto i = static_cast<decltype(y_.size())>(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<decltype(data_.size())>(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;
 
     }

include/boost/math/interpolators/septic_hermite.hpp

@@ -57,5 +57,27 @@ private:
     std::shared_ptr<detail::cardinal_septic_hermite_detail<RandomAccessContainer>> impl_;
 };
 
+
+template<class RandomAccessContainer>
+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<detail::cardinal_septic_hermite_detail_aos<RandomAccessContainer>>(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<detail::cardinal_septic_hermite_detail_aos<RandomAccessContainer>> impl_;
+};
+
 }
 #endif 

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<typename Real>
 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<std::array<Real, 4>> 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<typename Real>
 void test_linear()
 {
-    std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9};
+    std::vector<Real> x{0,1,2,3,4,5,6,7,8,9};
     std::vector<Real> y = x;
     std::vector<Real> dydx(x.size(), 1);
     std::vector<Real> 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<std::array<Real, 4>> 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<typename Real>
@@ -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<std::array<Real, 4>> 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<typename Real>
 void test_cubic()
 {