diff --git a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
index 224f7aab5..1df68f69e 100644
--- a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
+++ b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp
@@ -160,5 +160,225 @@ private:
     RandomAccessContainer d2ydx2_;
     RandomAccessContainer d3ydx3_;
 };
+
+template<class RandomAccessContainer>
+class cardinal_septic_hermite_detail {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    cardinal_septic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
+    : y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)}, x0_{x0}, dx_{dx}
+    {
+        if (y_.size() != dydx_.size()) {
+            throw std::domain_error("Numbers of derivatives must = number of ordinates.");
+        }
+        if (y_.size() != d2ydx2_.size()) {
+            throw std::domain_error("Number of second derivatives must equal number of ordinates.");
+        }
+        if (y_.size() != d3ydx3_.size()) {
+            throw std::domain_error("Number of third derivatives must equal number of ordinates.");
+        }
+
+        if (y_.size() < 2) {
+            throw std::domain_error("At least 2 abscissas are required.");
+        }
+    }
+
+    inline Real operator()(Real x) const {
+        Real xf = x0_ + (y_.size()-1)*dx_;
+        if  (x < x0_ || x > xf) {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::digits10+3);
+            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
+                << x0_ << ", " << xf << "]";
+            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 == xf)
+        {
+            return y_.back();
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const {
+        using std::floor;
+        auto i = static_cast<decltype(y_.size())>(floor((x-x0_)/dx_));
+        Real xi = x0_ + i*dx_;
+        Real t = (x - xi)/dx_;
+
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        // Velocity:
+        Real v0 = dydx_[i];
+        Real v1 = dydx_[i+1];
+        // Acceleration:
+        Real a0 = d2ydx2_[i];
+        Real a1 = d2ydx2_[i+1];
+        // Jerk:
+        Real j0 = d3ydx3_[i];
+        Real j1 = d3ydx3_[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 t5 = t4*t;
+        Real t6 = t5*t;
+        Real t7 = t6*t;
+        Real dx2 = dx_*dx_;
+        Real dx3 = dx2*dx_;
+
+        Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
+        Real z4 = -s;
+        Real z0 = s + 1;
+        Real z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t;
+        Real z2 = 2*t7 - 15*t6/2 + 10*t5 - 5*t4 + t2/2;
+        Real z3 = t7/6 - 2*t6/3 + t5 - 2*t4/3 + t3/6;
+
+        Real z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4;
+        Real z6 = -2*t7 + 13*t6/2 - 7*t5 + 5*t4/2;
+        Real z7 = t7/6 - t6/2 + t5/2 - t4/6;
+
+        return z0*y0 + z1*dx_*v0 + z2*dx2*a0 + z3*dx3*j0 + z4*y1 + z5*dx_*v1 + z6*dx2*a1 + z7*dx3*j1;
+    }
+
+    inline Real prime(Real x) const {
+        Real xf = x0_ + (y_.size()-1)*dx_;
+        if  (x < x0_ || x > xf)
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::digits10+3);
+            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
+                << x0_ << ", " << xf << "]";
+            throw std::domain_error(oss.str());
+        }
+        if (x == xf)
+        {
+            return dydx_.back();
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const {
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+
+private:
+    RandomAccessContainer y_;
+    RandomAccessContainer dydx_;
+    RandomAccessContainer d2ydx2_;
+    RandomAccessContainer d3ydx3_;
+    Real x0_;
+    Real dx_;
+};
+
+
+template<class RandomAccessContainer>
+class cardinal_septic_hermite_detail_aos {
+public:
+    using Point = typename RandomAccessContainer::value_type;
+    using Real = typename Point::value_type;
+    cardinal_septic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
+    : data_{std::move(dat)}, x0_{x0}, dx_{dx}
+    {
+        if (data_.size() < 2) {
+            throw std::domain_error("At least 2 abscissas are required.");
+        }
+        if (data_[0].size() != 4) {
+            throw std::domain_error("There must be 4 data items per struct.");
+        }
+    }
+
+    inline Real operator()(Real x) const {
+        Real xf = x0_ + (data_.size()-1)*dx_;
+        if  (x < x0_ || x > xf) {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::digits10+3);
+            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
+                << x0_ << ", " << xf << "]";
+            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 == xf)
+        {
+            return data_.back()[0];
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const {
+        using std::floor;
+        auto i = static_cast<decltype(data_.size())>(floor((x-x0_)/dx_));
+        Real xi = x0_ + i*dx_;
+        Real t = (x - xi)/dx_;
+
+        Real y0 = data_[i][0];
+        Real y1 = data_[i+1][0];
+        // Velocity:
+        Real v0 = data_[i][1];
+        Real v1 = 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;
+        Real t5 = t4*t;
+        Real t6 = t5*t;
+        Real t7 = t6*t;
+        Real dx2 = dx_*dx_;
+        Real dx3 = dx2*dx_;
+
+        Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
+        Real z4 = -s;
+        Real z0 = s + 1;
+        Real z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t;
+        Real z2 = 2*t7 - 15*t6/2 + 10*t5 - 5*t4 + t2/2;
+        Real z3 = t7/6 - 2*t6/3 + t5 - 2*t4/3 + t3/6;
+
+        Real z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4;
+        Real z6 = -2*t7 + 13*t6/2 - 7*t5 + 5*t4/2;
+        Real z7 = t7/6 - t6/2 + t5/2 - t4/6;
+
+        return z0*y0 + z1*dx_*v0 + z2*dx2*a0 + z3*dx3*j0 + z4*y1 + z5*dx_*v1 + z6*dx2*a1 + z7*dx3*j1;
+    }
+
+    inline Real prime(Real x) const {
+        Real xf = x0_ + (data_.size()-1)*dx_;
+        if  (x < x0_ || x > xf)
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::digits10+3);
+            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
+                << x0_ << ", " << xf << "]";
+            throw std::domain_error(oss.str());
+        }
+        if (x == xf)
+        {
+            return data_.back()[1];
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const {
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+
+private:
+    RandomAccessContainer data_;
+    Real x0_;
+    Real dx_;
+};
+
+
 }
-#endif 
\ No newline at end of file
+#endif
diff --git a/include/boost/math/interpolators/septic_hermite.hpp b/include/boost/math/interpolators/septic_hermite.hpp
index c544d1c63..54da49a02 100644
--- a/include/boost/math/interpolators/septic_hermite.hpp
+++ b/include/boost/math/interpolators/septic_hermite.hpp
@@ -34,5 +34,28 @@ public:
 private:
     std::shared_ptr<detail::septic_hermite_detail<RandomAccessContainer>> impl_;
 };
+
+template<class RandomAccessContainer>
+class cardinal_septic_hermite {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    cardinal_septic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx,
+                    RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
+     : impl_(std::make_shared<detail::cardinal_septic_hermite_detail<RandomAccessContainer>>(
+     std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), 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<RandomAccessContainer>> impl_;
+};
+
 }
 #endif 
\ No newline at end of file
diff --git a/include/boost/math/special_functions/daubechies_scaling.hpp b/include/boost/math/special_functions/daubechies_scaling.hpp
index 5416eb46f..43037e9f3 100644
--- a/include/boost/math/special_functions/daubechies_scaling.hpp
+++ b/include/boost/math/special_functions/daubechies_scaling.hpp
@@ -13,7 +13,9 @@
 #include <thread>
 #include <future>
 #include <iostream>
+#if BOOST_HAS_FLOAT128
 #include <boost/multiprecision/float128.hpp>
+#endif
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
 #include <boost/math/filters/daubechies.hpp>
@@ -89,7 +91,7 @@ public:
 
     matched_holder(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements) : y_{std::move(y)}, dydx_{std::move(dydx)}
     {
-        h_ = Real(1)/(1<< grid_refinements);
+        inv_h_ = 1<< grid_refinements;
     }
     
     inline Real operator()(Real x) const {
@@ -97,25 +99,27 @@ public:
         using std::floor;
         using std::sqrt;
         using std::pow;
-        if (x <= 0 || x >= 3) {
-            return 0;
-        }
         // This is the exact Holder exponent, but it's pessimistic almost everywhere!
         // It's only exactly right at dyadic rationals.
+        // Some compilers do not evaluate this at compile time:
         Real const alpha = 2 - log(1+sqrt(Real(3)))/log(Real(2));
+        Real const onemalpham1 = 1/(1-alpha);
+        //const Real alpha = 0.5500156865235042;
         // So we're gonna make the graph dip a little harder; this will capture more of the self-similar behavior:
         //Real constexpr const alpha = Real(3)/Real(10);
-        int64_t i = static_cast<int64_t>(floor(x/h_));
-        Real t = (x- i*h_)/h_;
+        Real s = x*inv_h_;
+        Real ii = floor(s);
+        int64_t i = static_cast<int64_t>(ii);
+        Real t = s - ii;
         Real v = y_[i];
-        Real dphi = dydx_[i+1]*h_;
-        v += (dphi - alpha*(y_[i+1] - y_[i]))*t/(1-alpha);
-        v += (-dphi + y_[i+1] - y_[i])*pow(t, alpha)/(1-alpha);
+        Real dphi = dydx_[i+1]/inv_h_;
+        Real diff = y_[i+1] - y_[i];
+        v += ((dphi - alpha*diff)*t + (-dphi + diff)*pow(t, alpha) )*onemalpham1;
         return v;
     }
 
 private:
-    Real h_;
+    Real inv_h_;
     RandomAccessContainer y_;
     RandomAccessContainer dydx_;
 };
@@ -155,7 +159,9 @@ public:
     daubechies_scaling(int grid_refinements = -1)
     {
         static_assert(p <= 15, "Scaling functions only implements up to p = 15");
+        #if BOOST_HAS_FLOAT128
         using boost::multiprecision::float128;
+        #endif
         if (grid_refinements < 0)
         {
             if (std::is_same_v<Real, float>)
@@ -275,13 +281,6 @@ public:
         auto y = t0.get();
         auto dydx = t1.get();
 
-        // Storing the vector of x's is unnecessary; it's only because I haven't implemented an equispaced cubic Hermite interpolator:
-        std::vector<Real> x(y.size());
-        Real h = Real(2*p-1)/(x.size()-1);
-        for (size_t i = 0; i < x.size(); ++i) {
-            x[i] = i*h;
-        }
-
         if constexpr (p==2) {
             m_mh = std::make_shared<detail::matched_holder<std::vector<Real>>>(std::move(y), std::move(dydx), grid_refinements);
         }
@@ -297,9 +296,9 @@ public:
             m_qh = std::make_shared<interpolators::detail::cardinal_quintic_hermite_detail<std::vector<Real>>>(std::move(y), std::move(dydx), std::move(d2ydx2), Real(0), dx);
         }
         if constexpr (p >= 10) {
-            m_sh = std::make_shared<interpolators::detail::septic_hermite_detail<std::vector<Real>>>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
+            Real dx = Real(1)/(1 << grid_refinements);
+            m_sh = std::make_shared<interpolators::detail::cardinal_septic_hermite_detail<std::vector<Real>>>(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), Real(0), dx);
         }
-
      }
 
 
@@ -320,7 +319,7 @@ public:
             return m_qh->unchecked_evaluation(x);
         }
         if constexpr (p >= 10) {
-            return m_sh->operator()(x);
+            return m_sh->unchecked_evaluation(x);
         }
     }
 
@@ -339,7 +338,7 @@ public:
             return m_qh->unchecked_prime(x);
         }
         if constexpr (p >= 10) {
-            return m_sh->prime(x);
+            return m_sh->unchecked_prime(x);
         }
     }
 
@@ -358,7 +357,7 @@ private:
     // Need this for p = 6,7,8,9:
     std::shared_ptr<interpolators::detail::cardinal_quintic_hermite_detail<std::vector<Real>>> m_qh;
     // Need this for p >= 10:
-    std::shared_ptr<interpolators::detail::septic_hermite_detail<std::vector<Real>>> m_sh;
+    std::shared_ptr<interpolators::detail::cardinal_septic_hermite_detail<std::vector<Real>>> m_sh;
 
     /*Real constant_interpolation(Real x) const {
         if (x <= 0 || x >= 2*p-1) {
diff --git a/test/septic_hermite_test.cpp b/test/septic_hermite_test.cpp
new file mode 100644
index 000000000..ed9772328
--- /dev/null
+++ b/test/septic_hermite_test.cpp
@@ -0,0 +1,288 @@
+/*
+ * 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 <numeric>
+#include <utility>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/math/interpolators/septic_hermite.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+using boost::multiprecision::float128;
+#endif
+
+
+using boost::math::interpolators::septic_hermite;
+using boost::math::interpolators::cardinal_septic_hermite;
+
+template<typename Real>
+void test_constant()
+{
+
+    std::vector<Real> x{0,1,2,3, 9, 22, 81};
+    std::vector<Real> y(x.size());
+    std::vector<Real> dydx(x.size(), 0);
+    std::vector<Real> d2ydx2(x.size(), 0);
+    std::vector<Real> d3ydx3(x.size(), 0);
+    for (auto & t : y) {
+        t = 7;
+    }
+
+    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) {
+        CHECK_ULP_CLOSE(Real(7), sh(t), 24);
+        CHECK_ULP_CLOSE(Real(0), sh.prime(t), 24);
+    }
+
+    Real x0 = 0;
+    Real dx = 1;
+    y.resize(128, 7);
+    dydx.resize(128, 0);
+    d2ydx2.resize(128, 0);
+    d3ydx3.resize(128, 0);
+    auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
+    for (Real t = x0; t <= 127; t += 0.25) {
+        CHECK_ULP_CLOSE(Real(7), csh(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> y = x;
+    std::vector<Real> dydx(x.size(), 1);
+    std::vector<Real> d2ydx2(x.size(), 0);
+    std::vector<Real> d3ydx3(x.size(), 0);
+
+    auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
+
+    for (Real t = 0; t <= 9; t += 0.25) {
+        CHECK_ULP_CLOSE(Real(t), sh(t), 2);
+        //CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2);
+    }
+
+    boost::random::mt19937 rng;
+    boost::random::uniform_real_distribution<Real> dis(0.5,1);
+    x.resize(512);
+    x[0] = dis(rng);
+    Real xmin = x[0];
+    for (size_t i = 1; i < x.size(); ++i) {
+        x[i] = x[i-1] + dis(rng);
+    }
+    Real xmax = x.back();
+
+    y = x;
+    dydx.resize(x.size(), 1);
+    d2ydx2.resize(x.size(), 0);
+    d3ydx3.resize(x.size(), 0);
+
+    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) {
+        CHECK_ULP_CLOSE(t, sh(t), 15);
+    }
+}
+
+template<typename Real>
+void test_quadratic()
+{
+    std::vector<Real> x{0,1,2,3,4,5,6,7,8,9};
+    std::vector<Real> y(x.size());
+    for (size_t i = 0; i < y.size(); ++i)
+    {
+        y[i] = x[i]*x[i]/2;
+    }
+
+    std::vector<Real> dydx(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        dydx[i] = x[i];
+    }
+
+    std::vector<Real> d2ydx2(x.size(), 1);
+    std::vector<Real> d3ydx3(x.size(), 0);
+
+    auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
+
+    for (Real t = 0; t <= 9; t += 0.0078125) {
+        CHECK_ULP_CLOSE(t*t/2, sh(t), 100);
+        //CHECK_ULP_CLOSE(t, qh.prime(t), 7);
+    }
+
+    boost::random::mt19937 rng;
+    boost::random::uniform_real_distribution<Real> dis(0.5,1);
+    x.resize(8);
+    x[0] = dis(rng);
+    Real xmin = x[0];
+    for (size_t i = 1; i < x.size(); ++i) {
+        x[i] = x[i-1] + dis(rng);
+    }
+    Real xmax = x.back();
+
+    y.resize(x.size());
+    for (size_t i = 0; i < y.size(); ++i)
+    {
+        y[i] = x[i]*x[i]/2;
+    }
+
+    dydx.resize(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        dydx[i] = x[i];
+    }
+
+    d2ydx2.resize(x.size(), 1);
+    d3ydx3.resize(x.size(), 0); 
+
+    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) {
+        CHECK_ULP_CLOSE(t*t/2, sh(t), 24);
+        //CHECK_ULP_CLOSE(t, qh.prime(t), 36);
+    }
+}
+
+template<typename Real>
+void test_cubic()
+{
+
+    std::vector<Real> x{0,1,2,3,4,5,6,7};
+    Real xmax = x.back();
+    std::vector<Real> y(x.size());
+    for (size_t i = 0; i < y.size(); ++i)
+    {
+        y[i] = x[i]*x[i]*x[i];
+    }
+
+    std::vector<Real> dydx(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        dydx[i] = 3*x[i]*x[i];
+    }
+
+    std::vector<Real> d2ydx2(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        d2ydx2[i] = 6*x[i];
+    }
+    std::vector<Real> d3ydx3(x.size(), 6);
+
+    auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
+
+    for (Real t = 0; t <= xmax; t += 0.0078125) {
+        CHECK_ULP_CLOSE(t*t*t, sh(t), 151);
+    }
+}
+
+template<typename Real>
+void test_quartic()
+{
+
+    std::vector<Real> x{0,1,2,3,4,5,6,7,8,9};
+    Real xmax = x.back();
+    std::vector<Real> y(x.size());
+    for (size_t i = 0; i < y.size(); ++i)
+    {
+        y[i] = x[i]*x[i]*x[i]*x[i];
+    }
+
+    std::vector<Real> dydx(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        dydx[i] = 4*x[i]*x[i]*x[i];
+    }
+
+    std::vector<Real> d2ydx2(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        d2ydx2[i] = 12*x[i]*x[i];
+    }
+
+    std::vector<Real> d3ydx3(x.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        d3ydx3[i] = 24*x[i];
+    }
+
+    auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
+
+    for (Real t = 1; t <= xmax; t += 0.0078125) {
+        CHECK_ULP_CLOSE(t*t*t*t, sh(t), 117);
+    }
+}
+
+
+template<typename Real>
+void test_interpolation_condition()
+{
+    for (size_t n = 4; n < 50; ++n) {
+        std::vector<Real> x(n);
+        std::vector<Real> y(n);
+        std::vector<Real> dydx(n);
+        std::vector<Real> d2ydx2(n);
+        std::vector<Real> d3ydx3(n);
+        boost::random::mt19937 rd; 
+        boost::random::uniform_real_distribution<Real> 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);
+            d3ydx3[i] = dis(rd);
+        }
+
+        auto x_copy = x;
+        auto y_copy = y;
+        auto dydx_copy = dydx;
+        auto d2ydx2_copy = d2ydx2;
+        auto d3ydx3_copy = d3ydx3;
+        auto s = septic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy), std::move(d3ydx3_copy));
+
+        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<float>();
+    test_linear<float>();
+    test_quadratic<float>();
+    test_cubic<float>();
+    test_quartic<float>();
+    test_interpolation_condition<float>();
+
+    test_constant<double>();
+    test_linear<double>();
+    test_quadratic<double>();
+    test_cubic<double>();
+    test_quartic<double>();
+    test_interpolation_condition<double>();
+
+    test_constant<long double>();
+    test_linear<long double>();
+    test_quadratic<long double>();
+    test_cubic<long double>();
+    test_quartic<long double>();
+    test_interpolation_condition<long double>();
+
+#ifdef BOOST_HAS_FLOAT128
+    test_constant<float128>();
+    test_linear<float128>();
+    test_quadratic<float128>();
+    test_cubic<float128>();
+    test_quartic<float128>();
+    test_interpolation_condition<float128>();
+#endif
+
+    return boost::math::test::report_errors();
+}