diff --git a/include/boost/math/interpolators/cubic_b_spline.hpp b/include/boost/math/interpolators/cubic_b_spline.hpp
index 73ac1d013..91a435aee 100644
--- a/include/boost/math/interpolators/cubic_b_spline.hpp
+++ b/include/boost/math/interpolators/cubic_b_spline.hpp
@@ -45,6 +45,8 @@ public:
 
     Real prime(Real x) const;
 
+    Real double_prime(Real x) const;
+
 private:
     std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
 };
@@ -74,5 +76,12 @@ Real cubic_b_spline<Real>::prime(Real x) const
     return m_imp->prime(x);
 }
 
+template<class Real>
+Real cubic_b_spline<Real>::double_prime(Real x) const
+{
+    return m_imp->double_prime(x);
+}
+
+
 }}
 #endif
diff --git a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
index de013fabb..1347cbc30 100644
--- a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
+++ b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
@@ -32,6 +32,8 @@ public:
 
     Real prime(Real x) const;
 
+    Real double_prime(Real x) const;
+
 private:
     std::vector<Real> m_beta;
     Real m_h_inv;
@@ -79,6 +81,25 @@ Real b3_spline_prime(Real x)
     return (Real) 0;
 }
 
+template<class Real>
+Real b3_spline_double_prime(Real x)
+{
+    if (x < 0)
+    {
+        return b3_spline_double_prime(-x);
+    }
+
+    if (x < 1)
+    {
+        return (3*boost::math::constants::half<Real>()*x - 2) + x*(3*boost::math::constants::half<Real>());
+    }
+    if (x < 2)
+    {
+        return (2 - x);
+    }
+    return (Real) 0;
+}
+
 
 template <class Real>
 template <class BidiIterator>
@@ -283,5 +304,27 @@ Real cubic_b_spline_imp<Real>::prime(Real x) const
     return z*m_h_inv;
 }
 
+template<class Real>
+Real cubic_b_spline_imp<Real>::double_prime(Real x) const
+{
+    Real z = 0;
+    Real t = m_h_inv*(x - m_a) + 1;
+
+    using std::max;
+    using std::min;
+    using std::ceil;
+    using std::floor;
+
+    size_t k_min = (size_t) (max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
+    size_t k_max = (size_t) (min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)));
+
+    for (size_t k = k_min; k <= k_max; ++k)
+    {
+        z += m_beta[k]*b3_spline_double_prime(t - k);
+    }
+    return z*m_h_inv*m_h_inv;
+}
+
+
 }}}
 #endif
diff --git a/test/test_cubic_b_spline.cpp b/test/test_cubic_b_spline.cpp
index 8655cabf2..8f4e66c65 100644
--- a/test/test_cubic_b_spline.cpp
+++ b/test/test_cubic_b_spline.cpp
@@ -30,6 +30,9 @@ void test_b3_spline()
     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(-2.5), (Real) 0);
     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(2.5), (Real) 0);
     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(-2.5), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime<Real>(2.5), (Real) 0);
+    BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime<Real>(-2.5), (Real) 0);
+
 
     // On the boundary of support:
     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(2), (Real) 0);
@@ -52,6 +55,7 @@ void test_b3_spline()
         Real arg = i*0.01;
         BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>(arg), boost::math::detail::b3_spline<Real>(arg), eps);
         BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime<Real>(-arg), -boost::math::detail::b3_spline_prime<Real>(arg), eps);
+        BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_double_prime<Real>(-arg), boost::math::detail::b3_spline_double_prime<Real>(arg), eps);
     }
 
 }
@@ -109,6 +113,9 @@ void test_constant_function()
         BOOST_CHECK_CLOSE(y, constant, 10*std::numeric_limits<Real>::epsilon());
         Real y_prime = spline.prime(i*step + a + 0.002);
         BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits<Real>::epsilon());
+        Real y_double_prime = spline.double_prime(i*step + a + 0.002);
+        BOOST_CHECK_SMALL(y_double_prime, 5000*std::numeric_limits<Real>::epsilon());
+
     }
 
     // Test that correctly specified left and right-derivatives work properly: