grabbed latest spline libs from BOOST, these are only in dev branch as of 4.1.2020

boost/math/interpolators/barycentric_rational.hpp

-/*
- *  Copyright Nick Thompson, 2017
- *  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)
- *
- *  Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an
- *  interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time.
- *  The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far.
- *  The measure of this stability is the "local mesh ratio", which can be queried from the routine.
- *  Pictorially, the following t_i spacing is bad (has a high local mesh ratio)
- *  ||             |      | |                           |
- *  and this t_i spacing is good (has a low local mesh ratio)
- *  |   |      |    |     |    |        |    |  |    |
- *
- *
- *  If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order.
- *  A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly.
- *
- *  References:
- *  Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation." Numerische Mathematik 107.2 (2007): 315-331.
- *  Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
- */
-
-#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
-#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
-
-#include <memory>
-#include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
-
-namespace boost{ namespace math{
-
-template<class Real>
-class barycentric_rational
-{
-public:
-    barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
-
-    barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
-
-    template <class InputIterator1, class InputIterator2>
-    barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type* = 0);
-
-    Real operator()(Real x) const;
-
-    Real prime(Real x) const;
-
-    std::vector<Real>&& return_x()
-    {
-        return m_imp->return_x();
-    }
-
-    std::vector<Real>&& return_y()
-    {
-        return m_imp->return_y();
-    }
-
-private:
-    std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
-};
-
-template <class Real>
-barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
- m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
-{
-    return;
-}
-
-template <class Real>
-barycentric_rational<Real>::barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order):
- m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(std::move(x), std::move(y), approximation_order))
-{
-    return;
-}
-
-
-template <class Real>
-template <class InputIterator1, class InputIterator2>
-barycentric_rational<Real>::barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type*)
- : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
-{
-}
-
-template<class Real>
-Real barycentric_rational<Real>::operator()(Real x) const
-{
-    return m_imp->operator()(x);
-}
-
-template<class Real>
-Real barycentric_rational<Real>::prime(Real x) const
-{
-    return m_imp->prime(x);
-}
-
-
-}}
-#endif
+/*
+ *  Copyright Nick Thompson, 2017
+ *  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)
+ *
+ *  Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an
+ *  interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time.
+ *  The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far.
+ *  The measure of this stability is the "local mesh ratio", which can be queried from the routine.
+ *  Pictorially, the following t_i spacing is bad (has a high local mesh ratio)
+ *  ||             |      | |                           |
+ *  and this t_i spacing is good (has a low local mesh ratio)
+ *  |   |      |    |     |    |        |    |  |    |
+ *
+ *
+ *  If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order.
+ *  A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly.
+ *
+ *  References:
+ *  Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation."
+*      Numerische Mathematik 107.2 (2007): 315-331.
+ *  Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
+
+#include <memory>
+#include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
+
+namespace boost{ namespace math{
+
+template<class Real>
+class barycentric_rational
+{
+public:
+    barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
+
+    barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
+
+    template <class InputIterator1, class InputIterator2>
+    barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type* = 0);
+
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    std::vector<Real>&& return_x()
+    {
+        return m_imp->return_x();
+    }
+
+    std::vector<Real>&& return_y()
+    {
+        return m_imp->return_y();
+    }
+
+private:
+    std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
+};
+
+template <class Real>
+barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
+ m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
+{
+    return;
+}
+
+template <class Real>
+barycentric_rational<Real>::barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order):
+ m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(std::move(x), std::move(y), approximation_order))
+{
+    return;
+}
+
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational<Real>::barycentric_rational(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order, typename boost::disable_if_c<boost::is_integral<InputIterator2>::value>::type*)
+ : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
+{
+}
+
+template<class Real>
+Real barycentric_rational<Real>::operator()(Real x) const
+{
+    return m_imp->operator()(x);
+}
+
+template<class Real>
+Real barycentric_rational<Real>::prime(Real x) const
+{
+    return m_imp->prime(x);
+}
+
+
+}}
+#endif

boost/math/interpolators/cardinal_cubic_b_spline.hpp

+// Copyright Nick Thompson, 2017
+// 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)
+
+// This implements the compactly supported cubic b spline algorithm described in
+// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
+// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
+
+// Let f be the function we are trying to interpolate, and s be the interpolating spline.
+// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
+// The order of accuracy depends on the regularity of the f, however, assuming f is
+// four-times continuously differentiable, the error is of O(h^4).
+// In addition, we can differentiate the spline and obtain a good interpolant for f'.
+// The main restriction of this method is that the samples of f must be evenly spaced.
+// Look for barycentric rational interpolation for non-evenly sampled data.
+// Properties:
+// - s(x_j) = f(x_j)
+// - All cubic polynomials interpolated exactly
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_HPP
+
+#include <boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp>
+
+namespace boost{ namespace math{ namespace interpolators {
+
+template <class Real>
+class cardinal_cubic_b_spline
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+    template <class BidiIterator>
+    cardinal_cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                   Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                   Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+    cardinal_cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+    cardinal_cubic_b_spline() = default;
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    Real double_prime(Real x) const;
+
+private:
+    std::shared_ptr<detail::cardinal_cubic_b_spline_imp<Real>> m_imp;
+};
+
+template<class Real>
+cardinal_cubic_b_spline<Real>::cardinal_cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+                                     Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cardinal_cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template <class Real>
+template <class BidiIterator>
+cardinal_cubic_b_spline<Real>::cardinal_cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+   Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cardinal_cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template<class Real>
+Real cardinal_cubic_b_spline<Real>::operator()(Real x) const
+{
+    return m_imp->operator()(x);
+}
+
+template<class Real>
+Real cardinal_cubic_b_spline<Real>::prime(Real x) const
+{
+    return m_imp->prime(x);
+}
+
+template<class Real>
+Real cardinal_cubic_b_spline<Real>::double_prime(Real x) const
+{
+    return m_imp->double_prime(x);
+}
+
+
+}}}
+#endif

boost/math/interpolators/cardinal_quadratic_b_spline.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp>
+
+
+namespace boost{ namespace math{ namespace interpolators {
+
+template <class Real>
+class cardinal_quadratic_b_spline
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // y[0] = y(a), y[n - 1] = y(b), step_size = (b - a)/(n -1).
+    cardinal_quadratic_b_spline(const Real* const y,
+                                size_t n,
+                                Real t0 /* initial time, left endpoint */,
+                                Real h  /*spacing, stepsize*/,
+                                Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                                Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
+     : impl_(std::make_shared<detail::cardinal_quadratic_b_spline_detail<Real>>(y, n, t0, h, left_endpoint_derivative, right_endpoint_derivative))
+    {}
+
+    // Oh the bizarre error messages if we template this on a RandomAccessContainer:
+    cardinal_quadratic_b_spline(std::vector<Real> const & y,
+                                Real t0 /* initial time, left endpoint */,
+                                Real h  /*spacing, stepsize*/,
+                                Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                                Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
+     : impl_(std::make_shared<detail::cardinal_quadratic_b_spline_detail<Real>>(y.data(), y.size(), t0, h, left_endpoint_derivative, right_endpoint_derivative))
+    {}
+
+
+    Real operator()(Real t) const {
+        return impl_->operator()(t);
+    }
+
+    Real prime(Real t) const {
+       return impl_->prime(t);
+    }
+
+    Real t_max() const {
+        return impl_->t_max();
+    }
+
+private:
+    std::shared_ptr<detail::cardinal_quadratic_b_spline_detail<Real>> impl_;
+};
+
+}}}
+#endif

boost/math/interpolators/cardinal_quintic_b_spline.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_HPP
+#include <memory>
+#include <limits>
+#include <boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp>
+
+
+namespace boost{ namespace math{ namespace interpolators {
+
+template <class Real>
+class cardinal_quintic_b_spline
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // y[0] = y(a), y[n - 1] = y(b), step_size = (b - a)/(n -1).
+    cardinal_quintic_b_spline(const Real* const y,
+                                size_t n,
+                                Real t0 /* initial time, left endpoint */,
+                                Real h  /*spacing, stepsize*/,
+                                std::pair<Real, Real> left_endpoint_derivatives = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN()},
+                                std::pair<Real, Real> right_endpoint_derivatives = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN()})
+     : impl_(std::make_shared<detail::cardinal_quintic_b_spline_detail<Real>>(y, n, t0, h, left_endpoint_derivatives, right_endpoint_derivatives))
+    {}
+
+    // Oh the bizarre error messages if we template this on a RandomAccessContainer:
+    cardinal_quintic_b_spline(std::vector<Real> const & y,
+                                Real t0 /* initial time, left endpoint */,
+                                Real h  /*spacing, stepsize*/,
+                                std::pair<Real, Real> left_endpoint_derivatives = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN()},
+                                std::pair<Real, Real> right_endpoint_derivatives = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN()})
+     : impl_(std::make_shared<detail::cardinal_quintic_b_spline_detail<Real>>(y.data(), y.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives))
+    {}
+
+
+    Real operator()(Real t) const {
+        return impl_->operator()(t);
+    }
+
+    Real prime(Real t) const {
+       return impl_->prime(t);
+    }
+
+    Real double_prime(Real t) const {
+        return impl_->double_prime(t);
+    }
+
+    Real t_max() const {
+        return impl_->t_max();
+    }
+
+private:
+    std::shared_ptr<detail::cardinal_quintic_b_spline_detail<Real>> impl_;
+};
+
+}}}
+#endif

boost/math/interpolators/cardinal_trigonometric.hpp

+//  (C) Copyright Nick Thompson 2019.
+//  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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_TRIGONOMETRIC_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_TRIGONOMETRIC_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/cardinal_trigonometric_detail.hpp>
+
+namespace boost { namespace math { namespace interpolators {
+
+template<class RandomAccessContainer>
+class cardinal_trigonometric
+{
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    cardinal_trigonometric(RandomAccessContainer const & v, Real t0, Real h)
+    {
+        m_impl = std::make_shared<interpolators::detail::cardinal_trigonometric_detail<Real>>(v.data(), v.size(), t0, h);
+    }
+
+    Real operator()(Real t) const
+    {
+        return m_impl->operator()(t);
+    }
+
+    Real prime(Real t) const
+    {
+        return m_impl->prime(t);
+    }
+
+    Real double_prime(Real t) const
+    {
+        return m_impl->double_prime(t);
+    }
+
+    Real period() const
+    {
+        return m_impl->period();
+    }
+
+    Real integrate() const
+    {
+        return m_impl->integrate();
+    }
+
+    Real squared_l2() const
+    {
+        return m_impl->squared_l2();
+    }
+
+private:
+    std::shared_ptr<interpolators::detail::cardinal_trigonometric_detail<Real>> m_impl;
+};
+
+}}}
+#endif

boost/math/interpolators/catmull_rom.hpp

-// Copyright Nick Thompson, 2017
-// 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)
-
-// This computes the Catmull-Rom spline from a list of points.
-
-#ifndef BOOST_MATH_INTERPOLATORS_CATMULL_ROM
-#define BOOST_MATH_INTERPOLATORS_CATMULL_ROM
-
-#include <cmath>
-#include <vector>
-#include <algorithm>
-#include <iterator>
-
-namespace boost{ namespace math{
-
-    namespace detail
-    {
-        template<class Point>
-        typename Point::value_type alpha_distance(Point const & p1, Point const & p2, typename Point::value_type alpha)
-        {
-            using std::pow;
-            using std::size;
-            typename Point::value_type dsq = 0;
-            for (size_t i = 0; i < size(p1); ++i)
-            {
-                typename Point::value_type dx = p1[i] - p2[i];
-                dsq += dx*dx;
-            }
-            return pow(dsq, alpha/2);
-        }
-    }
-
-template <class Point>
-class catmull_rom
-{
-public:
-
-    catmull_rom(std::vector<Point>&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2);
-
-    catmull_rom(std::initializer_list<Point> l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(std::vector<Point>(l), closed, alpha) {}
-
-    typename Point::value_type max_parameter() const
-    {
-        return m_max_s;
-    }
-
-    typename Point::value_type parameter_at_point(size_t i) const
-    {
-        return m_s[i+1];
-    }
-
-    Point operator()(const typename Point::value_type s) const;
-
-    Point prime(const typename Point::value_type s) const;
-
-    std::vector<Point>&& get_points()
-    {
-        return std::move(m_pnts);
-    }
-
-private:
-    std::vector<Point> m_pnts;
-    std::vector<typename Point::value_type> m_s;
-    typename Point::value_type m_max_s;
-};
-
-template<class Point>
-catmull_rom<Point>::catmull_rom(std::vector<Point>&& points, bool closed, typename Point::value_type alpha) : m_pnts(std::move(points))
-{
-    size_t num_pnts = m_pnts.size();
-    //std::cout << "Number of points = " << num_pnts << "\n";
-    if (num_pnts < 4)
-    {
-        throw std::domain_error("The Catmull-Rom curve requires at least 4 points.");
-    }
-    if (alpha < 0 || alpha > 1)
-    {
-        throw std::domain_error("The parametrization alpha must be in the range [0,1].");
-    }
-
-    using std::abs;
-    m_s.resize(num_pnts+3);
-    m_pnts.resize(num_pnts+3);
-    //std::cout << "Number of points now = " << m_pnts.size() << "\n";
-
-    m_pnts[num_pnts+1] = m_pnts[0];
-    m_pnts[num_pnts+2] = m_pnts[1];
-
-    auto tmp = m_pnts[num_pnts-1];
-    for (int64_t i = num_pnts-1; i >= 0; --i)
-    {
-        m_pnts[i+1] = m_pnts[i];
-    }
-    m_pnts[0] = tmp;
-
-    m_s[0] = -detail::alpha_distance<Point>(m_pnts[0], m_pnts[1], alpha);
-    if (abs(m_s[0]) < std::numeric_limits<typename Point::value_type>::epsilon())
-    {
-        throw std::domain_error("The first and last point should not be the same.\n");
-    }
-    m_s[1] = 0;
-    for (size_t i = 2; i < m_s.size(); ++i)
-    {
-        typename Point::value_type d = detail::alpha_distance<Point>(m_pnts[i], m_pnts[i-1], alpha);
-        if (abs(d) < std::numeric_limits<typename Point::value_type>::epsilon())
-        {
-            throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n");
-        }
-        m_s[i] = m_s[i-1] + d;
-    }
-    if(closed)
-    {
-        m_max_s = m_s[num_pnts+1];
-    }
-    else
-    {
-        m_max_s = m_s[num_pnts];
-    }
-}
-
-
-template<class Point>
-Point catmull_rom<Point>::operator()(const typename Point::value_type s) const
-{
-    using std::size;
-    if (s < 0 || s > m_max_s)
-    {
-        throw std::domain_error("Parameter outside bounds.");
-    }
-    auto it = std::upper_bound(m_s.begin(), m_s.end(), s);
-    //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]:
-    size_t i = std::distance(m_s.begin(), it - 1);
-
-    // Only denom21 is used twice:
-    typename Point::value_type denom21 = 1/(m_s[i+1] - m_s[i]);
-    typename Point::value_type s0s = m_s[i-1] - s;
-    typename Point::value_type s1s = m_s[i] - s;
-    typename Point::value_type s2s = m_s[i+1] - s;
-    typename Point::value_type s3s = m_s[i+2] - s;
-
-    Point A1_or_A3;
-    typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]);
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A1_or_A3[j] = denom*(s1s*m_pnts[i-1][j] - s0s*m_pnts[i][j]);
-    }
-
-    Point A2_or_B2;
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A2_or_B2[j] = denom21*(s2s*m_pnts[i][j] - s1s*m_pnts[i+1][j]);
-    }
-
-    Point B1_or_C;
-    denom = 1/(m_s[i+1] - m_s[i-1]);
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B1_or_C[j] = denom*(s2s*A1_or_A3[j] - s0s*A2_or_B2[j]);
-    }
-
-    denom = 1/(m_s[i+2] - m_s[i+1]);
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A1_or_A3[j] = denom*(s3s*m_pnts[i+1][j] - s2s*m_pnts[i+2][j]);
-    }
-
-    Point B2;
-    denom = 1/(m_s[i+2] - m_s[i]);
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B2[j] = denom*(s3s*A2_or_B2[j] - s1s*A1_or_A3[j]);
-    }
-
-    for(size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B1_or_C[j] = denom21*(s2s*B1_or_C[j] - s1s*B2[j]);
-    }
-
-    return B1_or_C;
-}
-
-template<class Point>
-Point catmull_rom<Point>::prime(const typename Point::value_type s) const
-{
-    using std::size;
-    // https://math.stackexchange.com/questions/843595/how-can-i-calculate-the-derivative-of-a-catmull-rom-spline-with-nonuniform-param
-    // http://denkovacs.com/2016/02/catmull-rom-spline-derivatives/
-    if (s < 0 || s > m_max_s)
-    {
-        throw std::domain_error("Parameter outside bounds.\n");
-    }
-    auto it = std::upper_bound(m_s.begin(), m_s.end(), s);
-    //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]:
-    size_t i = std::distance(m_s.begin(), it - 1);
-    Point A1;
-    typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]);
-    typename Point::value_type k1 = (m_s[i]-s)*denom;
-    typename Point::value_type k2 = (s - m_s[i-1])*denom;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j];
-    }
-
-    Point A1p;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A1p[j] = denom*(m_pnts[i][j] - m_pnts[i-1][j]);
-    }
-
-    Point A2;
-    denom = 1/(m_s[i+1] - m_s[i]);
-    k1 = (m_s[i+1]-s)*denom;
-    k2 = (s - m_s[i])*denom;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j];
-    }
-
-    Point A2p;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A2p[j] = denom*(m_pnts[i+1][j] - m_pnts[i][j]);
-    }
-
-
-    Point B1;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B1[j] = k1*A1[j] + k2*A2[j];
-    }
-
-    Point A3;
-    denom = 1/(m_s[i+2] - m_s[i+1]);
-    k1 = (m_s[i+2]-s)*denom;
-    k2 = (s - m_s[i+1])*denom;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j];
-    }
-
-    Point A3p;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        A3p[j] = denom*(m_pnts[i+2][j] - m_pnts[i+1][j]);
-    }
-
-    Point B2;
-    denom = 1/(m_s[i+2] - m_s[i]);
-    k1 = (m_s[i+2]-s)*denom;
-    k2 = (s - m_s[i])*denom;
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B2[j] = k1*A2[j] + k2*A3[j];
-    }
-
-    Point B1p;
-    denom = 1/(m_s[i+1] - m_s[i-1]);
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B1p[j] = denom*(A2[j] - A1[j] + (m_s[i+1]- s)*A1p[j] + (s-m_s[i-1])*A2p[j]);
-    }
-
-    Point B2p;
-    denom = 1/(m_s[i+2] - m_s[i]);
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        B2p[j] = denom*(A3[j] - A2[j] + (m_s[i+2] - s)*A2p[j] + (s - m_s[i])*A3p[j]);
-    }
-
-    Point Cp;
-    denom = 1/(m_s[i+1] - m_s[i]);
-    for (size_t j = 0; j < size(m_pnts[0]); ++j)
-    {
-        Cp[j] = denom*(B2[j] - B1[j] + (m_s[i+1] - s)*B1p[j] + (s - m_s[i])*B2p[j]);
-    }
-    return Cp;
-}
-
-
-}}
-#endif
+// Copyright Nick Thompson, 2017
+// 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)
+
+// This computes the Catmull-Rom spline from a list of points.
+
+#ifndef BOOST_MATH_INTERPOLATORS_CATMULL_ROM
+#define BOOST_MATH_INTERPOLATORS_CATMULL_ROM
+
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <iterator>
+#include <stdexcept>
+#include <boost/config.hpp>
+
+namespace std_workaround {
+
+#if defined(__cpp_lib_nonmember_container_access) || (defined(BOOST_MSVC) && (BOOST_MSVC >= 1900))
+   using std::size;
+#else
+   template <class C>
+   inline BOOST_CONSTEXPR std::size_t size(const C& c)
+   {
+      return c.size();
+   }
+   template <class T, std::size_t N>
+   inline BOOST_CONSTEXPR std::size_t size(const T(&array)[N]) BOOST_NOEXCEPT
+   {
+      return N;
+   }
+#endif
+}
+
+namespace boost{ namespace math{
+
+    namespace detail
+    {
+        template<class Point>
+        typename Point::value_type alpha_distance(Point const & p1, Point const & p2, typename Point::value_type alpha)
+        {
+            using std::pow;
+            using std_workaround::size;
+            typename Point::value_type dsq = 0;
+            for (size_t i = 0; i < size(p1); ++i)
+            {
+                typename Point::value_type dx = p1[i] - p2[i];
+                dsq += dx*dx;
+            }
+            return pow(dsq, alpha/2);
+        }
+    }
+
+template <class Point, class RandomAccessContainer = std::vector<Point> >
+class catmull_rom
+{
+   typedef typename Point::value_type value_type;
+public:
+
+    catmull_rom(RandomAccessContainer&& points, bool closed = false, value_type alpha = (value_type) 1/ (value_type) 2);
+
+    catmull_rom(std::initializer_list<Point> l, bool closed = false, value_type alpha = (value_type) 1/ (value_type) 2) : catmull_rom<Point, RandomAccessContainer>(RandomAccessContainer(l), closed, alpha) {}
+
+    value_type max_parameter() const
+    {
+        return m_max_s;
+    }
+
+    value_type parameter_at_point(size_t i) const
+    {
+        return m_s[i+1];
+    }
+
+    Point operator()(const value_type s) const;
+
+    Point prime(const value_type s) const;
+
+    RandomAccessContainer&& get_points()
+    {
+        return std::move(m_pnts);
+    }
+
+private:
+    RandomAccessContainer m_pnts;
+    std::vector<value_type> m_s;
+    value_type m_max_s;
+};
+
+template<class Point, class RandomAccessContainer >
+catmull_rom<Point, RandomAccessContainer>::catmull_rom(RandomAccessContainer&& points, bool closed, typename Point::value_type alpha) : m_pnts(std::move(points))
+{
+    std::size_t num_pnts = m_pnts.size();
+    //std::cout << "Number of points = " << num_pnts << "\n";
+    if (num_pnts < 4)
+    {
+        throw std::domain_error("The Catmull-Rom curve requires at least 4 points.");
+    }
+    if (alpha < 0 || alpha > 1)
+    {
+        throw std::domain_error("The parametrization alpha must be in the range [0,1].");
+    }
+
+    using std::abs;
+    m_s.resize(num_pnts+3);
+    m_pnts.resize(num_pnts+3);
+    //std::cout << "Number of points now = " << m_pnts.size() << "\n";
+
+    m_pnts[num_pnts+1] = m_pnts[0];
+    m_pnts[num_pnts+2] = m_pnts[1];
+
+    auto tmp = m_pnts[num_pnts-1];
+    for (std::ptrdiff_t i = num_pnts-1; i >= 0; --i)
+    {
+        m_pnts[i+1] = m_pnts[i];
+    }
+    m_pnts[0] = tmp;
+
+    m_s[0] = -detail::alpha_distance<Point>(m_pnts[0], m_pnts[1], alpha);
+    if (abs(m_s[0]) < std::numeric_limits<typename Point::value_type>::epsilon())
+    {
+        throw std::domain_error("The first and last point should not be the same.\n");
+    }
+    m_s[1] = 0;
+    for (size_t i = 2; i < m_s.size(); ++i)
+    {
+        typename Point::value_type d = detail::alpha_distance<Point>(m_pnts[i], m_pnts[i-1], alpha);
+        if (abs(d) < std::numeric_limits<typename Point::value_type>::epsilon())
+        {
+            throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n");
+        }
+        m_s[i] = m_s[i-1] + d;
+    }
+    if(closed)
+    {
+        m_max_s = m_s[num_pnts+1];
+    }
+    else
+    {
+        m_max_s = m_s[num_pnts];
+    }
+}
+
+
+template<class Point, class RandomAccessContainer >
+Point catmull_rom<Point, RandomAccessContainer>::operator()(const typename Point::value_type s) const
+{
+    using std_workaround::size;
+    if (s < 0 || s > m_max_s)
+    {
+        throw std::domain_error("Parameter outside bounds.");
+    }
+    auto it = std::upper_bound(m_s.begin(), m_s.end(), s);
+    //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]:
+    size_t i = std::distance(m_s.begin(), it - 1);
+
+    // Only denom21 is used twice:
+    typename Point::value_type denom21 = 1/(m_s[i+1] - m_s[i]);
+    typename Point::value_type s0s = m_s[i-1] - s;
+    typename Point::value_type s1s = m_s[i] - s;
+    typename Point::value_type s2s = m_s[i+1] - s;
+    typename Point::value_type s3s = m_s[i+2] - s;
+
+    Point A1_or_A3;
+    typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]);
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A1_or_A3[j] = denom*(s1s*m_pnts[i-1][j] - s0s*m_pnts[i][j]);
+    }
+
+    Point A2_or_B2;
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A2_or_B2[j] = denom21*(s2s*m_pnts[i][j] - s1s*m_pnts[i+1][j]);
+    }
+
+    Point B1_or_C;
+    denom = 1/(m_s[i+1] - m_s[i-1]);
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B1_or_C[j] = denom*(s2s*A1_or_A3[j] - s0s*A2_or_B2[j]);
+    }
+
+    denom = 1/(m_s[i+2] - m_s[i+1]);
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A1_or_A3[j] = denom*(s3s*m_pnts[i+1][j] - s2s*m_pnts[i+2][j]);
+    }
+
+    Point B2;
+    denom = 1/(m_s[i+2] - m_s[i]);
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B2[j] = denom*(s3s*A2_or_B2[j] - s1s*A1_or_A3[j]);
+    }
+
+    for(size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B1_or_C[j] = denom21*(s2s*B1_or_C[j] - s1s*B2[j]);
+    }
+
+    return B1_or_C;
+}
+
+template<class Point, class RandomAccessContainer >
+Point catmull_rom<Point, RandomAccessContainer>::prime(const typename Point::value_type s) const
+{
+    using std_workaround::size;
+    // https://math.stackexchange.com/questions/843595/how-can-i-calculate-the-derivative-of-a-catmull-rom-spline-with-nonuniform-param
+    // http://denkovacs.com/2016/02/catmull-rom-spline-derivatives/
+    if (s < 0 || s > m_max_s)
+    {
+        throw std::domain_error("Parameter outside bounds.\n");
+    }
+    auto it = std::upper_bound(m_s.begin(), m_s.end(), s);
+    //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]:
+    size_t i = std::distance(m_s.begin(), it - 1);
+    Point A1;
+    typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]);
+    typename Point::value_type k1 = (m_s[i]-s)*denom;
+    typename Point::value_type k2 = (s - m_s[i-1])*denom;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j];
+    }
+
+    Point A1p;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A1p[j] = denom*(m_pnts[i][j] - m_pnts[i-1][j]);
+    }
+
+    Point A2;
+    denom = 1/(m_s[i+1] - m_s[i]);
+    k1 = (m_s[i+1]-s)*denom;
+    k2 = (s - m_s[i])*denom;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j];
+    }
+
+    Point A2p;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A2p[j] = denom*(m_pnts[i+1][j] - m_pnts[i][j]);
+    }
+
+
+    Point B1;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B1[j] = k1*A1[j] + k2*A2[j];
+    }
+
+    Point A3;
+    denom = 1/(m_s[i+2] - m_s[i+1]);
+    k1 = (m_s[i+2]-s)*denom;
+    k2 = (s - m_s[i+1])*denom;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j];
+    }
+
+    Point A3p;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        A3p[j] = denom*(m_pnts[i+2][j] - m_pnts[i+1][j]);
+    }
+
+    Point B2;
+    denom = 1/(m_s[i+2] - m_s[i]);
+    k1 = (m_s[i+2]-s)*denom;
+    k2 = (s - m_s[i])*denom;
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B2[j] = k1*A2[j] + k2*A3[j];
+    }
+
+    Point B1p;
+    denom = 1/(m_s[i+1] - m_s[i-1]);
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B1p[j] = denom*(A2[j] - A1[j] + (m_s[i+1]- s)*A1p[j] + (s-m_s[i-1])*A2p[j]);
+    }
+
+    Point B2p;
+    denom = 1/(m_s[i+2] - m_s[i]);
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        B2p[j] = denom*(A3[j] - A2[j] + (m_s[i+2] - s)*A2p[j] + (s - m_s[i])*A3p[j]);
+    }
+
+    Point Cp;
+    denom = 1/(m_s[i+1] - m_s[i]);
+    for (size_t j = 0; j < size(m_pnts[0]); ++j)
+    {
+        Cp[j] = denom*(B2[j] - B1[j] + (m_s[i+1] - s)*B1p[j] + (s - m_s[i])*B2p[j]);
+    }
+    return Cp;
+}
+
+
+}}
+#endif

boost/math/interpolators/cubic_b_spline.hpp

-// Copyright Nick Thompson, 2017
-// 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)
-
-// This implements the compactly supported cubic b spline algorithm described in
-// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
-// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
-
-// Let f be the function we are trying to interpolate, and s be the interpolating spline.
-// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
-// The order of accuracy depends on the regularity of the f, however, assuming f is
-// four-times continuously differentiable, the error is of O(h^4).
-// In addition, we can differentiate the spline and obtain a good interpolant for f'.
-// The main restriction of this method is that the samples of f must be evenly spaced.
-// Look for barycentric rational interpolation for non-evenly sampled data.
-// Properties:
-// - s(x_j) = f(x_j)
-// - All cubic polynomials interpolated exactly
-
-#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
-#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
-
-#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
-
-namespace boost{ namespace math{
-
-template <class Real>
-class cubic_b_spline
-{
-public:
-    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
-    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
-    template <class BidiIterator>
-    cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
-                   Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
-                   Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
-    cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
-       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
-       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
-
-    cubic_b_spline() = default;
-    Real operator()(Real x) const;
-
-    Real prime(Real x) const;
-
-private:
-    std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
-};
-
-template<class Real>
-cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
-                                     Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
-{
-}
-
-template <class Real>
-template <class BidiIterator>
-cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
-   Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
-{
-}
-
-template<class Real>
-Real cubic_b_spline<Real>::operator()(Real x) const
-{
-    return m_imp->operator()(x);
-}
-
-template<class Real>
-Real cubic_b_spline<Real>::prime(Real x) const
-{
-    return m_imp->prime(x);
-}
-
-}}
-#endif
+// Copyright Nick Thompson, 2017
+// 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)
+
+// This implements the compactly supported cubic b spline algorithm described in
+// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
+// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
+
+// Let f be the function we are trying to interpolate, and s be the interpolating spline.
+// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
+// The order of accuracy depends on the regularity of the f, however, assuming f is
+// four-times continuously differentiable, the error is of O(h^4).
+// In addition, we can differentiate the spline and obtain a good interpolant for f'.
+// The main restriction of this method is that the samples of f must be evenly spaced.
+// Look for barycentric rational interpolation for non-evenly sampled data.
+// Properties:
+// - s(x_j) = f(x_j)
+// - All cubic polynomials interpolated exactly
+
+#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
+
+#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
+#include <boost/config/header_deprecated.hpp>
+
+BOOST_HEADER_DEPRECATED("<boost/math/interpolators/cardinal_cubic_b_spline.hpp>");
+
+namespace boost{ namespace math{
+
+template <class Real>
+class cubic_b_spline
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+    template <class BidiIterator>
+    cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                   Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                   Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+    cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+    cubic_b_spline() = default;
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    Real double_prime(Real x) const;
+
+private:
+    std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
+};
+
+template<class Real>
+cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
+                                     Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+   Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
+{
+}
+
+template<class Real>
+Real cubic_b_spline<Real>::operator()(Real x) const
+{
+    return m_imp->operator()(x);
+}
+
+template<class Real>
+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

boost/math/interpolators/cubic_hermite.hpp

+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_HERMITE_HPP
+#define BOOST_MATH_INTERPOLATORS_CUBIC_HERMITE_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template<class RandomAccessContainer>
+class cubic_hermite {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    cubic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx) 
+    : impl_(std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(dydx)))
+    {}
+
+    inline Real operator()(Real x) const {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const {
+        return impl_->prime(x);
+    }
+
+    friend std::ostream& operator<<(std::ostream & os, const cubic_hermite & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    void push_back(Real x, Real y, Real dydx)
+    {
+        impl_->push_back(x, y, dydx);
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+private:
+    std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+template<class RandomAccessContainer>
+class cardinal_cubic_hermite {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    cardinal_cubic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, Real x0, Real dx) 
+    : impl_(std::make_shared<detail::cardinal_cubic_hermite_detail<RandomAccessContainer>>(std::move(y), std::move(dydx), x0, dx))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    friend std::ostream& operator<<(std::ostream & os, const cardinal_cubic_hermite & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+private:
+    std::shared_ptr<detail::cardinal_cubic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+
+template<class RandomAccessContainer>
+class cardinal_cubic_hermite_aos {
+public:
+    using Point = typename RandomAccessContainer::value_type;
+    using Real = typename Point::value_type;
+
+    cardinal_cubic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx) 
+    : impl_(std::make_shared<detail::cardinal_cubic_hermite_detail_aos<RandomAccessContainer>>(std::move(data), x0, dx))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    friend std::ostream& operator<<(std::ostream & os, const cardinal_cubic_hermite_aos & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+private:
+    std::shared_ptr<detail::cardinal_cubic_hermite_detail_aos<RandomAccessContainer>> impl_;
+};
+
+
+}
+#endif
\ No newline at end of file

boost/math/interpolators/detail/barycentric_rational_detail.hpp

-/*
- *  Copyright Nick Thompson, 2017
- *  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)
- */
-
-#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
-#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
-
-#include <vector>
-#include <utility> // for std::move
-#include <algorithm> // for std::is_sorted
-#include <boost/lexical_cast.hpp>
-#include <boost/math/special_functions/fpclassify.hpp>
-#include <boost/core/demangle.hpp>
-#include <boost/assert.hpp>
-
-namespace boost{ namespace math{ namespace detail{
-
-template<class Real>
-class barycentric_rational_imp
-{
-public:
-    template <class InputIterator1, class InputIterator2>
-    barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
-
-    barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
-
-    Real operator()(Real x) const;
-
-    Real prime(Real x) const;
-
-    // The barycentric weights are not really that interesting; except to the unit tests!
-    Real weight(size_t i) const { return m_w[i]; }
-
-    std::vector<Real>&& return_x()
-    {
-        return std::move(m_x);
-    }
-
-    std::vector<Real>&& return_y()
-    {
-        return std::move(m_y);
-    }
-
-private:
-
-    void calculate_weights(size_t approximation_order);
-
-    std::vector<Real> m_x;
-    std::vector<Real> m_y;
-    std::vector<Real> m_w;
-};
-
-template <class Real>
-template <class InputIterator1, class InputIterator2>
-barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
-{
-    std::ptrdiff_t n = std::distance(start_x, end_x);
-
-    if (approximation_order >= (std::size_t)n)
-    {
-        throw std::domain_error("Approximation order must be < data length.");
-    }
-
-    // Big sad memcpy.
-    m_x.resize(n);
-    m_y.resize(n);
-    for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
-    {
-        // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
-        if(boost::math::isnan(*start_x))
-        {
-            std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
-            throw std::domain_error(msg);
-        }
-
-        if(boost::math::isnan(*start_y))
-        {
-           std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
-           throw std::domain_error(msg);
-        }
-
-        m_x[i] = *start_x;
-        m_y[i] = *start_y;
-    }
-    calculate_weights(approximation_order);
-}
-
-template <class Real>
-barycentric_rational_imp<Real>::barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y))
-{
-    BOOST_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates.");
-    BOOST_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length.");
-    BOOST_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1].");
-    calculate_weights(approximation_order);
-}
-
-template<class Real>
-void barycentric_rational_imp<Real>::calculate_weights(size_t approximation_order)
-{
-    using std::abs;
-    int64_t n = m_x.size();
-    m_w.resize(n, 0);
-    for(int64_t k = 0; k < n; ++k)
-    {
-        int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
-        int64_t i_max = k;
-        if (k >= n - (std::ptrdiff_t)approximation_order)
-        {
-            i_max = n - approximation_order - 1;
-        }
-
-        for(int64_t i = i_min; i <= i_max; ++i)
-        {
-            Real inv_product = 1;
-            int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
-            for(int64_t j = i; j <= j_max; ++j)
-            {
-                if (j == k)
-                {
-                    continue;
-                }
-
-                Real diff = m_x[k] - m_x[j];
-                using std::numeric_limits;
-                if (abs(diff) < (numeric_limits<Real>::min)())
-                {
-                   std::string msg = std::string("Spacing between  x[")
-                      + boost::lexical_cast<std::string>(k) + std::string("] and x[")
-                      + boost::lexical_cast<std::string>(i) + std::string("] is ")
-                      + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ")
-                      + boost::core::demangle(typeid(Real).name());
-                    throw std::logic_error(msg);
-                }
-                inv_product *= diff;
-            }
-            if (i % 2 == 0)
-            {
-                m_w[k] += 1/inv_product;
-            }
-            else
-            {
-                m_w[k] -= 1/inv_product;
-            }
-        }
-    }
-}
-
-
-template<class Real>
-Real barycentric_rational_imp<Real>::operator()(Real x) const
-{
-    Real numerator = 0;
-    Real denominator = 0;
-    for(size_t i = 0; i < m_x.size(); ++i)
-    {
-        // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
-        // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
-        // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
-        if (x == m_x[i])
-        {
-            return m_y[i];
-        }
-        Real t = m_w[i]/(x - m_x[i]);
-        numerator += t*m_y[i];
-        denominator += t;
-    }
-    return numerator/denominator;
-}
-
-/*
- * A formula for computing the derivative of the barycentric representation is given in
- * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
- * Mathematics of Computation, v47, number 175, 1986.
- * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
- * and reviewed in
- * Recent developments in barycentric rational interpolation
- * Jean–Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
- *
- * Is it possible to complete this in one pass through the data?
- */
-
-template<class Real>
-Real barycentric_rational_imp<Real>::prime(Real x) const
-{
-    Real rx = this->operator()(x);
-    Real numerator = 0;
-    Real denominator = 0;
-    for(size_t i = 0; i < m_x.size(); ++i)
-    {
-        if (x == m_x[i])
-        {
-            Real sum = 0;
-            for (size_t j = 0; j < m_x.size(); ++j)
-            {
-                if (j == i)
-                {
-                    continue;
-                }
-                sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]);
-            }
-            return -sum/m_w[i];
-        }
-        Real t = m_w[i]/(x - m_x[i]);
-        Real diff = (rx - m_y[i])/(x-m_x[i]);
-        numerator += t*diff;
-        denominator += t;
-    }
-
-    return numerator/denominator;
-}
-}}}
-#endif
+/*
+ *  Copyright Nick Thompson, 2017
+ *  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)
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
+
+#include <vector>
+#include <utility> // for std::move
+#include <algorithm> // for std::is_sorted
+#include <boost/lexical_cast.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/core/demangle.hpp>
+#include <boost/assert.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template<class Real>
+class barycentric_rational_imp
+{
+public:
+    template <class InputIterator1, class InputIterator2>
+    barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
+
+    barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
+
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    // The barycentric weights are not really that interesting; except to the unit tests!
+    Real weight(size_t i) const { return m_w[i]; }
+
+    std::vector<Real>&& return_x()
+    {
+        return std::move(m_x);
+    }
+
+    std::vector<Real>&& return_y()
+    {
+        return std::move(m_y);
+    }
+
+private:
+
+    void calculate_weights(size_t approximation_order);
+
+    std::vector<Real> m_x;
+    std::vector<Real> m_y;
+    std::vector<Real> m_w;
+};
+
+template <class Real>
+template <class InputIterator1, class InputIterator2>
+barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
+{
+    std::ptrdiff_t n = std::distance(start_x, end_x);
+
+    if (approximation_order >= (std::size_t)n)
+    {
+        throw std::domain_error("Approximation order must be < data length.");
+    }
+
+    // Big sad memcpy.
+    m_x.resize(n);
+    m_y.resize(n);
+    for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
+    {
+        // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
+        if(boost::math::isnan(*start_x))
+        {
+            std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+            throw std::domain_error(msg);
+        }
+
+        if(boost::math::isnan(*start_y))
+        {
+           std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
+           throw std::domain_error(msg);
+        }
+
+        m_x[i] = *start_x;
+        m_y[i] = *start_y;
+    }
+    calculate_weights(approximation_order);
+}
+
+template <class Real>
+barycentric_rational_imp<Real>::barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y))
+{
+    BOOST_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates.");
+    BOOST_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length.");
+    BOOST_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1].");
+    calculate_weights(approximation_order);
+}
+
+template<class Real>
+void barycentric_rational_imp<Real>::calculate_weights(size_t approximation_order)
+{
+    using std::abs;
+    int64_t n = m_x.size();
+    m_w.resize(n, 0);
+    for(int64_t k = 0; k < n; ++k)
+    {
+        int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
+        int64_t i_max = k;
+        if (k >= n - (std::ptrdiff_t)approximation_order)
+        {
+            i_max = n - approximation_order - 1;
+        }
+
+        for(int64_t i = i_min; i <= i_max; ++i)
+        {
+            Real inv_product = 1;
+            int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
+            for(int64_t j = i; j <= j_max; ++j)
+            {
+                if (j == k)
+                {
+                    continue;
+                }
+
+                Real diff = m_x[k] - m_x[j];
+                using std::numeric_limits;
+                if (abs(diff) < (numeric_limits<Real>::min)())
+                {
+                   std::string msg = std::string("Spacing between  x[")
+                      + boost::lexical_cast<std::string>(k) + std::string("] and x[")
+                      + boost::lexical_cast<std::string>(i) + std::string("] is ")
+                      + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ")
+                      + boost::core::demangle(typeid(Real).name());
+                    throw std::logic_error(msg);
+                }
+                inv_product *= diff;
+            }
+            if (i % 2 == 0)
+            {
+                m_w[k] += 1/inv_product;
+            }
+            else
+            {
+                m_w[k] -= 1/inv_product;
+            }
+        }
+    }
+}
+
+
+template<class Real>
+Real barycentric_rational_imp<Real>::operator()(Real x) const
+{
+    Real numerator = 0;
+    Real denominator = 0;
+    for(size_t i = 0; i < m_x.size(); ++i)
+    {
+        // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
+        // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
+        // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
+        if (x == m_x[i])
+        {
+            return m_y[i];
+        }
+        Real t = m_w[i]/(x - m_x[i]);
+        numerator += t*m_y[i];
+        denominator += t;
+    }
+    return numerator/denominator;
+}
+
+/*
+ * A formula for computing the derivative of the barycentric representation is given in
+ * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
+ * Mathematics of Computation, v47, number 175, 1986.
+ * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
+ * and reviewed in
+ * Recent developments in barycentric rational interpolation
+ * Jean-Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
+ *
+ * Is it possible to complete this in one pass through the data?
+ */
+
+template<class Real>
+Real barycentric_rational_imp<Real>::prime(Real x) const
+{
+    Real rx = this->operator()(x);
+    Real numerator = 0;
+    Real denominator = 0;
+    for(size_t i = 0; i < m_x.size(); ++i)
+    {
+        if (x == m_x[i])
+        {
+            Real sum = 0;
+            for (size_t j = 0; j < m_x.size(); ++j)
+            {
+                if (j == i)
+                {
+                    continue;
+                }
+                sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]);
+            }
+            return -sum/m_w[i];
+        }
+        Real t = m_w[i]/(x - m_x[i]);
+        Real diff = (rx - m_y[i])/(x-m_x[i]);
+        numerator += t*diff;
+        denominator += t;
+    }
+
+    return numerator/denominator;
+}
+}}}
+#endif

boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp

+// Copyright Nick Thompson, 2017
+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP
+
+#include <limits>
+#include <cmath>
+#include <vector>
+#include <memory>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{ namespace math{ namespace interpolators{ namespace detail{
+
+
+template <class Real>
+class cardinal_cubic_b_spline_imp
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+    template <class BidiIterator>
+    cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    Real double_prime(Real x) const;
+
+private:
+    std::vector<Real> m_beta;
+    Real m_h_inv;
+    Real m_a;
+    Real m_avg;
+};
+
+
+
+template <class Real>
+Real b3_spline(Real x)
+{
+    using std::abs;
+    Real absx = abs(x);
+    if (absx < 1)
+    {
+        Real y = 2 - absx;
+        Real z = 1 - absx;
+        return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
+    }
+    if (absx < 2)
+    {
+        Real y = 2 - absx;
+        return boost::math::constants::sixth<Real>()*y*y*y;
+    }
+    return (Real) 0;
+}
+
+template<class Real>
+Real b3_spline_prime(Real x)
+{
+    if (x < 0)
+    {
+        return -b3_spline_prime(-x);
+    }
+
+    if (x < 1)
+    {
+        return x*(3*boost::math::constants::half<Real>()*x - 2);
+    }
+    if (x < 2)
+    {
+        return -boost::math::constants::half<Real>()*(2 - x)*(2 - 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*x - 2;
+    }
+    if (x < 2)
+    {
+        return (2 - x);
+    }
+    return (Real) 0;
+}
+
+
+template <class Real>
+template <class BidiIterator>
+cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                                             Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
+{
+    using boost::math::constants::third;
+
+    std::size_t length = end_p - f;
+
+    if (length < 5)
+    {
+        if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
+        {
+            throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
+        }
+        if (length < 3)
+        {
+            throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
+        }
+    }
+
+    if (boost::math::isnan(left_endpoint))
+    {
+        throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
+    }
+    if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
+    {
+        throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
+    }
+    if (step_size <= 0)
+    {
+        throw std::logic_error("The step size must be strictly > 0.\n");
+    }
+    // Storing the inverse of the stepsize does provide a measurable speedup.
+    // It's not huge, but nonetheless worthwhile.
+    m_h_inv = 1/step_size;
+
+    // Following Kress's notation, s'(a) = a1, s'(b) = b1
+    Real a1 = left_endpoint_derivative;
+    // See the finite-difference table on Wikipedia for reference on how
+    // to construct high-order estimates for one-sided derivatives:
+    // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
+    // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
+    if (boost::math::isnan(a1))
+    {
+        // For simple functions (linear, quadratic, so on)
+        // almost all the error comes from derivative estimation.
+        // This does pairwise summation which gives us another digit of accuracy over naive summation.
+        Real t0 = 4*(f[1] + third<Real>()*f[3]);
+        Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
+        a1 = m_h_inv*(t0 + t1);
+    }
+
+    Real b1 = right_endpoint_derivative;
+    if (boost::math::isnan(b1))
+    {
+        size_t n = length - 1;
+        Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
+        Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4  - 3*f[n - 2];
+
+        b1 = m_h_inv*(t0 + t1);
+    }
+
+    // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
+    // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
+    m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+    // Since the splines have compact support, they decay to zero very fast outside the endpoints.
+    // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
+    // boundary [a,b] without massive error.
+    // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
+    // This algorithm for computing the average is recommended in
+    // http://www.heikohoffmann.de/htmlthesis/node134.html
+    Real t = 1;
+    for (size_t i = 0; i < length; ++i)
+    {
+        if (boost::math::isnan(f[i]))
+        {
+            std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
+            throw std::logic_error(err);
+        }
+        m_avg += (f[i] - m_avg) / t;
+        t += 1;
+    }
+
+
+    // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
+    // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
+    // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
+    // See Kress, equations 8.41
+    // The the "tridiagonal" matrix is:
+    // 1  0 -1
+    // 1  4  1
+    //    1  4  1
+    //       1  4  1
+    //          ....
+    //          1  4  1
+    //          1  0 -1
+    // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
+    std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
+    std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+    rhs[0] = -2*step_size*a1;
+    rhs[rhs.size() - 1] = -2*step_size*b1;
+
+    super_diagonal[0] = 0;
+
+    for(size_t i = 1; i < rhs.size() - 1; ++i)
+    {
+        rhs[i] = 6*(f[i - 1] - m_avg);
+        super_diagonal[i] = 1;
+    }
+
+
+    // One step of row reduction on the first row to patch up the 5-diagonal problem:
+    // 1 0 -1 | r0
+    // 1 4 1  | r1
+    // mapsto:
+    // 1 0 -1 | r0
+    // 0 4 2  | r1 - r0
+    // mapsto
+    // 1 0 -1 | r0
+    // 0 1 1/2| (r1 - r0)/4
+    super_diagonal[1] = 0.5;
+    rhs[1] = (rhs[1] - rhs[0])/4;
+
+    // Now do a tridiagonal row reduction the standard way, until just before the last row:
+    for (size_t i = 2; i < rhs.size() - 1; ++i)
+    {
+        Real diagonal = 4 - super_diagonal[i - 1];
+        rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
+        super_diagonal[i] /= diagonal;
+    }
+
+    // Now the last row, which is in the form
+    // 1 sd[n-3] 0      | rhs[n-3]
+    // 0  1     sd[n-2] | rhs[n-2]
+    // 1  0     -1      | rhs[n-1]
+    Real final_subdiag = -super_diagonal[rhs.size() - 3];
+    rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
+    Real final_diag = -1/final_subdiag;
+    // Now we're here:
+    // 1 sd[n-3] 0         | rhs[n-3]
+    // 0  1     sd[n-2]    | rhs[n-2]
+    // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
+
+    final_diag = final_diag - super_diagonal[rhs.size() - 2];
+    rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
+
+
+    // Back substitutions:
+    m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
+    for(size_t i = rhs.size() - 2; i > 0; --i)
+    {
+        m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
+    }
+    m_beta[0] = m_beta[2] + rhs[0];
+}
+
+template<class Real>
+Real cardinal_cubic_b_spline_imp<Real>::operator()(Real x) const
+{
+    // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
+    // just the (at most 5) whose support overlaps the argument.
+    Real z = m_avg;
+    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) (max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0);
+
+    for (size_t k = k_min; k <= k_max; ++k)
+    {
+        z += m_beta[k]*b3_spline(t - k);
+    }
+
+    return z;
+}
+
+template<class Real>
+Real cardinal_cubic_b_spline_imp<Real>::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_prime(t - k);
+    }
+    return z*m_h_inv;
+}
+
+template<class Real>
+Real cardinal_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

boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
+#include <vector>
+#include <cmath>
+#include <stdexcept>
+
+namespace boost{ namespace math{ namespace interpolators{ namespace detail{
+
+template <class Real>
+Real b2_spline(Real x) {
+    using std::abs;
+    Real absx = abs(x);
+    if (absx < 1/Real(2))
+    {
+        Real y = absx - 1/Real(2);
+        Real z = absx + 1/Real(2);
+        return (2-y*y-z*z)/2;
+    }
+    if (absx < Real(3)/Real(2))
+    {
+        Real y = absx - Real(3)/Real(2);
+        return y*y/2;
+    }
+    return (Real) 0;
+}
+
+template <class Real>
+Real b2_spline_prime(Real x) {
+    if (x < 0) {
+        return -b2_spline_prime(-x);
+    }
+
+    if (x < 1/Real(2))
+    {
+        return -2*x;
+    }
+    if (x < Real(3)/Real(2))
+    {
+        return x - Real(3)/Real(2);
+    }
+    return (Real) 0;
+}
+
+
+template <class Real>
+class cardinal_quadratic_b_spline_detail
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
+
+    cardinal_quadratic_b_spline_detail(const Real* const y,
+                                size_t n,
+                                Real t0 /* initial time, left endpoint */,
+                                Real h  /*spacing, stepsize*/,
+                                Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                                Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
+    {
+        if (h <= 0) {
+            throw std::logic_error("Spacing must be > 0.");
+        }
+        m_inv_h = 1/h;
+        m_t0 = t0;
+
+        if (n < 3) {
+            throw std::logic_error("The interpolator requires at least 3 points.");
+        }
+
+        using std::isnan;
+        Real a;
+        if (isnan(left_endpoint_derivative)) {
+            // http://web.media.mit.edu/~crtaylor/calculator.html
+            a = -3*y[0] + 4*y[1] - y[2];
+        }
+        else {
+            a = 2*h*left_endpoint_derivative;
+        }
+
+        Real b;
+        if (isnan(right_endpoint_derivative)) {
+            b = 3*y[n-1] - 4*y[n-2] + y[n-3];
+        }
+        else {
+            b = 2*h*right_endpoint_derivative;
+        }
+
+        m_alpha.resize(n + 2);
+
+        // Begin row reduction:
+        std::vector<Real> rhs(n + 2, std::numeric_limits<Real>::quiet_NaN());
+        std::vector<Real> super_diagonal(n + 2, std::numeric_limits<Real>::quiet_NaN());
+
+        rhs[0] = -a;
+        rhs[rhs.size() - 1] = b;
+
+        super_diagonal[0] = 0;
+
+        for(size_t i = 1; i < rhs.size() - 1; ++i) {
+            rhs[i] = 8*y[i - 1];
+            super_diagonal[i] = 1;
+        }
+
+        // Patch up 5-diagonal problem:
+        rhs[1] = (rhs[1] - rhs[0])/6;
+        super_diagonal[1] = Real(1)/Real(3);
+        // First two rows are now:
+        // 1 0 -1 | -2hy0'
+        // 0 1 1/3| (8y0+2hy0')/6
+
+
+        // Start traditional tridiagonal row reduction:
+        for (size_t i = 2; i < rhs.size() - 1; ++i) {
+            Real diagonal = 6 - super_diagonal[i - 1];
+            rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
+            super_diagonal[i] /= diagonal;
+        }
+
+        //  1 sd[n-1] 0     | rhs[n-1]
+        //  0 1       sd[n] | rhs[n]
+        // -1 0       1     | rhs[n+1]
+
+        rhs[n+1] = rhs[n+1] + rhs[n-1];
+        Real bottom_subdiagonal = super_diagonal[n-1];
+
+        // We're here:
+        //  1 sd[n-1] 0     | rhs[n-1]
+        //  0 1       sd[n] | rhs[n]
+        //  0 bs      1     | rhs[n+1]
+
+        rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]);
+
+        m_alpha[n+1] = rhs[n+1];
+        for (size_t i = n; i > 0; --i) {
+            m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i];
+        }
+        m_alpha[0] = m_alpha[2] + rhs[0];
+    }
+
+    Real operator()(Real t) const {
+        if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
+            const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
+            throw std::domain_error(err_msg);
+        }
+        // Let k, gamma be defined via t = t0 + kh + gamma * h.
+        // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
+        // j_min = ceil((t-t0)/h - 1/2)
+        // j_max = floor(t-t0)/h + 5/2)
+        using std::floor;
+        using std::ceil;
+        Real x = (t-m_t0)*m_inv_h;
+        size_t j_min = ceil(x - Real(1)/Real(2));
+        size_t j_max = ceil(x + Real(5)/Real(2));
+        if (j_max >= m_alpha.size()) {
+            j_max = m_alpha.size() - 1;
+        }
+
+        Real y = 0;
+        x += 1;
+        for (size_t j = j_min; j <= j_max; ++j) {
+            y += m_alpha[j]*detail::b2_spline(x - j);
+        }
+        return y;
+    }
+
+    Real prime(Real t) const {
+        if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
+            const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
+            throw std::domain_error(err_msg);
+        }
+        // Let k, gamma be defined via t = t0 + kh + gamma * h.
+        // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
+        // j_min = ceil((t-t0)/h - 1/2)
+        // j_max = floor(t-t0)/h + 5/2)
+        using std::floor;
+        using std::ceil;
+        Real x = (t-m_t0)*m_inv_h;
+        size_t j_min = ceil(x - Real(1)/Real(2));
+        size_t j_max = ceil(x + Real(5)/Real(2));
+        if (j_max >= m_alpha.size()) {
+            j_max = m_alpha.size() - 1;
+        }
+
+        Real y = 0;
+        x += 1;
+        for (size_t j = j_min; j <= j_max; ++j) {
+            y += m_alpha[j]*detail::b2_spline_prime(x - j);
+        }
+        return y*m_inv_h;
+    }
+
+    Real t_max() const {
+        return m_t0 + (m_alpha.size()-3)/m_inv_h;
+    }
+
+private:
+    std::vector<Real> m_alpha;
+    Real m_inv_h;
+    Real m_t0;
+};
+
+}}}}
+#endif

boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
+#include <cmath>
+#include <vector>
+#include <utility>
+#include <boost/math/special_functions/cardinal_b_spline.hpp>
+
+namespace boost{ namespace math{ namespace interpolators{ namespace detail{
+
+
+template <class Real>
+class cardinal_quintic_b_spline_detail
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
+
+    cardinal_quintic_b_spline_detail(const Real* const y,
+                                     size_t n,
+                                     Real t0 /* initial time, left endpoint */,
+                                     Real h  /*spacing, stepsize*/,
+                                     std::pair<Real, Real> left_endpoint_derivatives,
+                                     std::pair<Real, Real> right_endpoint_derivatives)
+    {
+        static_assert(!std::is_integral<Real>::value, "The quintic B-spline interpolator only works with floating point types.");
+        if (h <= 0) {
+            throw std::logic_error("Spacing must be > 0.");
+        }
+        m_inv_h = 1/h;
+        m_t0 = t0;
+
+        if (n < 8) {
+            throw std::logic_error("The quintic B-spline interpolator requires at least 8 points.");
+        }
+
+        using std::isnan;
+        // This interpolator has error of order h^6, so the derivatives should be estimated with the same error.
+        // See: https://en.wikipedia.org/wiki/Finite_difference_coefficient
+        if (isnan(left_endpoint_derivatives.first)) {
+            Real tmp = -49*y[0]/20 + 6*y[1] - 15*y[2]/2 + 20*y[3]/3 - 15*y[4]/4 + 6*y[5]/5 - y[6]/6;
+            left_endpoint_derivatives.first = tmp/h;
+        }
+        if (isnan(right_endpoint_derivatives.first)) {
+            Real tmp = 49*y[n-1]/20 - 6*y[n-2] + 15*y[n-3]/2 - 20*y[n-4]/3 + 15*y[n-5]/4 - 6*y[n-6]/5 + y[n-7]/6;
+            right_endpoint_derivatives.first = tmp/h;
+        }
+        if(isnan(left_endpoint_derivatives.second)) {
+            Real tmp = 469*y[0]/90 - 223*y[1]/10 + 879*y[2]/20 - 949*y[3]/18 + 41*y[4] - 201*y[5]/10 + 1019*y[6]/180 - 7*y[7]/10;
+            left_endpoint_derivatives.second = tmp/(h*h);
+        }
+        if (isnan(right_endpoint_derivatives.second)) {
+            Real tmp = 469*y[n-1]/90 - 223*y[n-2]/10 + 879*y[n-3]/20 - 949*y[n-4]/18 + 41*y[n-5] - 201*y[n-6]/10 + 1019*y[n-7]/180 - 7*y[n-8]/10;
+            right_endpoint_derivatives.second = tmp/(h*h);
+        }
+
+        // This is really challenging my mental limits on by-hand row reduction.
+        // I debated bringing in a dependency on a sparse linear solver, but given that that would cause much agony for users I decided against it.
+
+        std::vector<Real> rhs(n+4);
+        rhs[0] = 20*y[0] - 12*h*left_endpoint_derivatives.first +  2*h*h*left_endpoint_derivatives.second;
+        rhs[1] = 60*y[0] - 12*h*left_endpoint_derivatives.first;
+        for (size_t i = 2; i < n + 2; ++i) {
+            rhs[i] = 120*y[i-2];
+        }
+        rhs[n+2] = 60*y[n-1] + 12*h*right_endpoint_derivatives.first;
+        rhs[n+3] = 20*y[n-1] + 12*h*right_endpoint_derivatives.first +  2*h*h*right_endpoint_derivatives.second;
+
+        std::vector<Real> diagonal(n+4, 66);
+        diagonal[0] = 1;
+        diagonal[1] = 18;
+        diagonal[n+2] = 18;
+        diagonal[n+3] = 1;
+
+        std::vector<Real> first_superdiagonal(n+4, 26);
+        first_superdiagonal[0] = 10;
+        first_superdiagonal[1] = 33;
+        first_superdiagonal[n+2] = 1;
+        // There is one less superdiagonal than diagonal; make sure that if we read it, it shows up as a bug:
+        first_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
+
+        std::vector<Real> second_superdiagonal(n+4, 1);
+        second_superdiagonal[0] = 9;
+        second_superdiagonal[1] = 8;
+        second_superdiagonal[n+2] = std::numeric_limits<Real>::quiet_NaN();
+        second_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
+
+        std::vector<Real> first_subdiagonal(n+4, 26);
+        first_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
+        first_subdiagonal[1] = 1;
+        first_subdiagonal[n+2] = 33;
+        first_subdiagonal[n+3] = 10;
+
+        std::vector<Real> second_subdiagonal(n+4, 1);
+        second_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
+        second_subdiagonal[1] = std::numeric_limits<Real>::quiet_NaN();
+        second_subdiagonal[n+2] = 8;
+        second_subdiagonal[n+3] = 9;
+
+
+        for (size_t i = 0; i < n+2; ++i) {
+            Real di = diagonal[i];
+            diagonal[i] = 1;
+            first_superdiagonal[i] /= di;
+            second_superdiagonal[i] /= di;
+            rhs[i] /= di;
+
+            // Eliminate first subdiagonal:
+            Real nfsub = -first_subdiagonal[i+1];
+            // Superfluous:
+            first_subdiagonal[i+1] /= nfsub;
+            // Not superfluous:
+            diagonal[i+1] /= nfsub;
+            first_superdiagonal[i+1] /= nfsub;
+            second_superdiagonal[i+1] /= nfsub;
+            rhs[i+1] /= nfsub;
+
+            diagonal[i+1] += first_superdiagonal[i];
+            first_superdiagonal[i+1] += second_superdiagonal[i];
+            rhs[i+1] += rhs[i];
+            // Superfluous, but clarifying:
+            first_subdiagonal[i+1] = 0;
+
+            // Eliminate second subdiagonal:
+            Real nssub = -second_subdiagonal[i+2];
+            first_subdiagonal[i+2] /= nssub;
+            diagonal[i+2] /= nssub;
+            first_superdiagonal[i+2] /= nssub;
+            second_superdiagonal[i+2] /= nssub;
+            rhs[i+2] /= nssub;
+
+            first_subdiagonal[i+2] += first_superdiagonal[i];
+            diagonal[i+2] += second_superdiagonal[i];
+            rhs[i+2] += rhs[i];
+            // Superfluous, but clarifying:
+            second_subdiagonal[i+2] = 0;
+        }
+
+        // Eliminate last subdiagonal:
+        Real dnp2 = diagonal[n+2];
+        diagonal[n+2] = 1;
+        first_superdiagonal[n+2] /= dnp2;
+        rhs[n+2] /= dnp2;
+        Real nfsubnp3 = -first_subdiagonal[n+3];
+        diagonal[n+3] /= nfsubnp3;
+        rhs[n+3] /= nfsubnp3;
+
+        diagonal[n+3] += first_superdiagonal[n+2];
+        rhs[n+3] += rhs[n+2];
+
+        m_alpha.resize(n + 4, std::numeric_limits<Real>::quiet_NaN());
+
+        m_alpha[n+3] = rhs[n+3]/diagonal[n+3];
+        m_alpha[n+2] = rhs[n+2] - first_superdiagonal[n+2]*m_alpha[n+3];
+        for (int64_t i = int64_t(n+1); i >= 0; --i) {
+            m_alpha[i] = rhs[i] - first_superdiagonal[i]*m_alpha[i+1] - second_superdiagonal[i]*m_alpha[i+2];
+        }
+
+    }
+
+    Real operator()(Real t) const {
+        using std::ceil;
+        using std::floor;
+        using boost::math::cardinal_b_spline;
+        // tf = t0 + (n-1)*h
+        // alpha.size() = n+4
+        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
+            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
+            throw std::domain_error(err_msg);
+        }
+        Real x = (t-m_t0)*m_inv_h;
+        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5.
+        // TODO: Zero pad m_alpha so that only the domain check is necessary.
+        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
+        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
+        Real s = 0;
+        for (int64_t j = j_min; j <= j_max; ++j) {
+            // TODO: Use Cox 1972 to generate all integer translates of B5 simultaneously.
+            s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2);
+        }
+        return s;
+    }
+
+    Real prime(Real t) const {
+        using std::ceil;
+        using std::floor;
+        using boost::math::cardinal_b_spline_prime;
+        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
+            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
+            throw std::domain_error(err_msg);
+        }
+        Real x = (t-m_t0)*m_inv_h;
+        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
+        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
+        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
+        Real s = 0;
+        for (int64_t j = j_min; j <= j_max; ++j) {
+            s += m_alpha[j]*cardinal_b_spline_prime<5, Real>(x - j + 2);
+        }
+        return s*m_inv_h;
+
+    }
+
+    Real double_prime(Real t) const {
+        using std::ceil;
+        using std::floor;
+        using boost::math::cardinal_b_spline_double_prime;
+        if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
+            const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
+            throw std::domain_error(err_msg);
+        }
+        Real x = (t-m_t0)*m_inv_h;
+        // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
+        int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1)));
+        int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
+        Real s = 0;
+        for (int64_t j = j_min; j <= j_max; ++j) {
+            s += m_alpha[j]*cardinal_b_spline_double_prime<5, Real>(x - j + 2);
+        }
+        return s*m_inv_h*m_inv_h;
+    }
+
+
+    Real t_max() const {
+        return m_t0 + (m_alpha.size()-5)/m_inv_h;
+    }
+
+private:
+    std::vector<Real> m_alpha;
+    Real m_inv_h;
+    Real m_t0;
+};
+
+}}}}
+#endif

boost/math/interpolators/detail/cardinal_trigonometric_detail.hpp

+//  (C) Copyright Nick Thompson 2019.
+//  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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_CARDINAL_TRIGONOMETRIC_HPP
+#define BOOST_MATH_INTERPOLATORS_DETAIL_CARDINAL_TRIGONOMETRIC_HPP
+#include <cmath>
+#include <stdexcept>
+#include <fftw3.h>
+#include <boost/math/constants/constants.hpp>
+
+#ifdef BOOST_HAS_FLOAT128
+#include <quadmath.h>
+#endif
+
+namespace boost { namespace math { namespace interpolators { namespace detail {
+
+template<typename Real>
+class cardinal_trigonometric_detail {
+public:
+  cardinal_trigonometric_detail(const Real* data, size_t length, Real t0, Real h)
+  {
+    m_data = data;
+    m_length = length;
+    m_t0 = t0;
+    m_h = h;
+    throw std::domain_error("Not implemented.");
+  }
+private:
+  size_t m_length;
+  Real m_t0;
+  Real m_h;
+  Real* m_data;
+};
+
+template<>
+class cardinal_trigonometric_detail<float> {
+public:
+  cardinal_trigonometric_detail(const float* data, size_t length, float t0, float h) : m_t0{t0}, m_h{h}
+  {
+    if (length == 0)
+    {
+      throw std::logic_error("At least one sample is required.");
+    }
+    if (h <= 0)
+    {
+      throw std::logic_error("The step size must be > 0");
+    }
+    // The period sadly must be stored, since the complex vector has length that cannot be used to recover the period:
+    m_T = m_h*length;
+    m_complex_vector_size = length/2 + 1;
+    m_gamma = fftwf_alloc_complex(m_complex_vector_size);
+    // The const_cast is legitimate: FFTW does not change the data as long as FFTW_ESTIMATE is provided.
+    fftwf_plan plan = fftwf_plan_dft_r2c_1d(length, const_cast<float*>(data), m_gamma, FFTW_ESTIMATE);
+    // FFTW says a null plan is impossible with the basic interface we are using, and I have no reason to doubt them.
+    // But it just feels weird not to check this:
+    if (!plan)
+    {
+      throw std::logic_error("A null fftw plan was created.");
+    }
+
+    fftwf_execute(plan);
+    fftwf_destroy_plan(plan);
+
+    float denom = length;
+    for (size_t k = 0; k < m_complex_vector_size; ++k)
+    {
+      m_gamma[k][0] /= denom;
+      m_gamma[k][1] /= denom;
+    }
+
+    if (length % 2 == 0)
+    {
+      m_gamma[m_complex_vector_size -1][0] /= 2;
+      // numerically, m_gamma[m_complex_vector_size -1][1] should be zero . . .
+      // I believe, but need to check, that FFTW guarantees that it is identically zero.
+    }
+  }
+
+  cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old)  = delete;
+
+  cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete;
+
+  cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete;
+
+  float operator()(float t) const
+  {
+    using std::sin;
+    using std::cos;
+    using boost::math::constants::two_pi;
+    using std::exp;
+    float s = m_gamma[0][0];
+    float x = two_pi<float>()*(t - m_t0)/m_T;
+    fftwf_complex z;
+    // boost::math::cos_pi with a redefinition of x? Not now . . .
+    z[0] = cos(x);
+    z[1] = sin(x);
+    fftwf_complex b{0, 0};
+    // u = b*z
+    fftwf_complex u;
+    for (size_t k = m_complex_vector_size - 1; k >= 1; --k) {
+      u[0] = b[0]*z[0] - b[1]*z[1];
+      u[1] = b[0]*z[1] + b[1]*z[0];
+      b[0] = m_gamma[k][0] + u[0];
+      b[1] = m_gamma[k][1] + u[1];
+    }
+
+    s += 2*(b[0]*z[0] - b[1]*z[1]);
+    return s;
+  }
+
+  float prime(float t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      float x = two_pi<float>()*(t - m_t0)/m_T;
+      fftwf_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftwf_complex b{0, 0};
+      // u = b*z
+      fftwf_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*m_gamma[k][0] + u[0];
+        b[1] = k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<float>()*(b[1]*z[0] + b[0]*z[1])/m_T;
+  }
+
+  float double_prime(float t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      float x = two_pi<float>()*(t - m_t0)/m_T;
+      fftwf_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftwf_complex b{0, 0};
+      // u = b*z
+      fftwf_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*k*m_gamma[k][0] + u[0];
+        b[1] = k*k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<float>()*two_pi<float>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T);
+  }
+
+  float period() const
+  {
+    return m_T;
+  }
+
+  float integrate() const
+  {
+    return m_T*m_gamma[0][0];
+  }
+
+  float squared_l2() const
+  {
+    float s = 0;
+    // Always add smallest to largest for accuracy.
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
+
+  ~cardinal_trigonometric_detail()
+  {
+    if (m_gamma)
+    {
+      fftwf_free(m_gamma);
+      m_gamma = nullptr;
+    }
+  }
+
+
+private:
+  float m_t0;
+  float m_h;
+  float m_T;
+  fftwf_complex* m_gamma;
+  size_t m_complex_vector_size;
+};
+
+
+template<>
+class cardinal_trigonometric_detail<double> {
+public:
+  cardinal_trigonometric_detail(const double* data, size_t length, double t0, double h) : m_t0{t0}, m_h{h}
+  {
+    if (length == 0)
+    {
+      throw std::logic_error("At least one sample is required.");
+    }
+    if (h <= 0)
+    {
+      throw std::logic_error("The step size must be > 0");
+    }
+    m_T = m_h*length;
+    m_complex_vector_size = length/2 + 1;
+    m_gamma = fftw_alloc_complex(m_complex_vector_size);
+    fftw_plan plan = fftw_plan_dft_r2c_1d(length, const_cast<double*>(data), m_gamma, FFTW_ESTIMATE);
+    if (!plan)
+    {
+      throw std::logic_error("A null fftw plan was created.");
+    }
+
+    fftw_execute(plan);
+    fftw_destroy_plan(plan);
+
+    double denom = length;
+    for (size_t k = 0; k < m_complex_vector_size; ++k)
+    {
+      m_gamma[k][0] /= denom;
+      m_gamma[k][1] /= denom;
+    }
+
+    if (length % 2 == 0)
+    {
+      m_gamma[m_complex_vector_size -1][0] /= 2;
+    }
+  }
+
+  cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old)  = delete;
+
+  cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete;
+
+  cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete;
+
+  double operator()(double t) const
+  {
+    using std::sin;
+    using std::cos;
+    using boost::math::constants::two_pi;
+    using std::exp;
+    double s = m_gamma[0][0];
+    double x = two_pi<double>()*(t - m_t0)/m_T;
+    fftw_complex z;
+    z[0] = cos(x);
+    z[1] = sin(x);
+    fftw_complex b{0, 0};
+    // u = b*z
+    fftw_complex u;
+    for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+    {
+      u[0] = b[0]*z[0] - b[1]*z[1];
+      u[1] = b[0]*z[1] + b[1]*z[0];
+      b[0] = m_gamma[k][0] + u[0];
+      b[1] = m_gamma[k][1] + u[1];
+    }
+
+    s += 2*(b[0]*z[0] - b[1]*z[1]);
+    return s;
+  }
+
+  double prime(double t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      double x = two_pi<double>()*(t - m_t0)/m_T;
+      fftw_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftw_complex b{0, 0};
+      // u = b*z
+      fftw_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*m_gamma[k][0] + u[0];
+        b[1] = k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<double>()*(b[1]*z[0] + b[0]*z[1])/m_T;
+  }
+
+  double double_prime(double t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      double x = two_pi<double>()*(t - m_t0)/m_T;
+      fftw_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftw_complex b{0, 0};
+      // u = b*z
+      fftw_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*k*m_gamma[k][0] + u[0];
+        b[1] = k*k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<double>()*two_pi<double>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T);
+  }
+
+  double period() const
+  {
+    return m_T;
+  }
+
+  double integrate() const
+  {
+    return m_T*m_gamma[0][0];
+  }
+
+  double squared_l2() const
+  {
+    double s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
+  ~cardinal_trigonometric_detail()
+  {
+    if (m_gamma)
+    {
+      fftw_free(m_gamma);
+      m_gamma = nullptr;
+    }
+  }
+
+private:
+  double m_t0;
+  double m_h;
+  double m_T;
+  fftw_complex* m_gamma;
+  size_t m_complex_vector_size;
+};
+
+
+template<>
+class cardinal_trigonometric_detail<long double> {
+public:
+  cardinal_trigonometric_detail(const long double* data, size_t length, long double t0, long double h) : m_t0{t0}, m_h{h}
+  {
+    if (length == 0)
+    {
+      throw std::logic_error("At least one sample is required.");
+    }
+    if (h <= 0)
+    {
+      throw std::logic_error("The step size must be > 0");
+    }
+    m_T = m_h*length;
+    m_complex_vector_size = length/2 + 1;
+    m_gamma = fftwl_alloc_complex(m_complex_vector_size);
+    fftwl_plan plan = fftwl_plan_dft_r2c_1d(length, const_cast<long double*>(data), m_gamma, FFTW_ESTIMATE);
+    if (!plan)
+    {
+      throw std::logic_error("A null fftw plan was created.");
+    }
+
+    fftwl_execute(plan);
+    fftwl_destroy_plan(plan);
+
+    long double denom = length;
+    for (size_t k = 0; k < m_complex_vector_size; ++k)
+    {
+      m_gamma[k][0] /= denom;
+      m_gamma[k][1] /= denom;
+    }
+
+    if (length % 2 == 0) {
+      m_gamma[m_complex_vector_size -1][0] /= 2;
+    }
+  }
+
+  cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old)  = delete;
+
+  cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete;
+
+  cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete;
+
+  long double operator()(long double t) const
+  {
+    using std::sin;
+    using std::cos;
+    using boost::math::constants::two_pi;
+    using std::exp;
+    long double s = m_gamma[0][0];
+    long double x = two_pi<long double>()*(t - m_t0)/m_T;
+    fftwl_complex z;
+    z[0] = cos(x);
+    z[1] = sin(x);
+    fftwl_complex b{0, 0};
+    fftwl_complex u;
+    for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+    {
+      u[0] = b[0]*z[0] - b[1]*z[1];
+      u[1] = b[0]*z[1] + b[1]*z[0];
+      b[0] = m_gamma[k][0] + u[0];
+      b[1] = m_gamma[k][1] + u[1];
+    }
+
+    s += 2*(b[0]*z[0] - b[1]*z[1]);
+    return s;
+  }
+
+  long double prime(long double t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      long double x = two_pi<long double>()*(t - m_t0)/m_T;
+      fftwl_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftwl_complex b{0, 0};
+      // u = b*z
+      fftwl_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*m_gamma[k][0] + u[0];
+        b[1] = k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<long double>()*(b[1]*z[0] + b[0]*z[1])/m_T;
+  }
+
+  long double double_prime(long double t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      long double x = two_pi<long double>()*(t - m_t0)/m_T;
+      fftwl_complex z;
+      z[0] = cos(x);
+      z[1] = sin(x);
+      fftwl_complex b{0, 0};
+      // u = b*z
+      fftwl_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*k*m_gamma[k][0] + u[0];
+        b[1] = k*k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<long double>()*two_pi<long double>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T);
+  }
+
+  long double period() const
+  {
+    return m_T;
+  }
+
+  long double integrate() const
+  {
+    return m_T*m_gamma[0][0];
+  }
+
+  long double squared_l2() const
+  {
+    long double s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
+  ~cardinal_trigonometric_detail()
+  {
+    if (m_gamma)
+    {
+      fftwl_free(m_gamma);
+      m_gamma = nullptr;
+    }
+  }
+
+private:
+  long double m_t0;
+  long double m_h;
+  long double m_T;
+  fftwl_complex* m_gamma;
+  size_t m_complex_vector_size;
+};
+
+#ifdef BOOST_HAS_FLOAT128
+template<>
+class cardinal_trigonometric_detail<__float128> {
+public:
+  cardinal_trigonometric_detail(const __float128* data, size_t length, __float128 t0, __float128 h) : m_t0{t0}, m_h{h}
+  {
+    if (length == 0)
+    {
+      throw std::logic_error("At least one sample is required.");
+    }
+    if (h <= 0)
+    {
+      throw std::logic_error("The step size must be > 0");
+    }
+    m_T = m_h*length;
+    m_complex_vector_size = length/2 + 1;
+    m_gamma = fftwq_alloc_complex(m_complex_vector_size);
+    fftwq_plan plan = fftwq_plan_dft_r2c_1d(length, reinterpret_cast<__float128*>(const_cast<__float128*>(data)), m_gamma, FFTW_ESTIMATE);
+    if (!plan)
+    {
+      throw std::logic_error("A null fftw plan was created.");
+    }
+
+    fftwq_execute(plan);
+    fftwq_destroy_plan(plan);
+
+    __float128 denom = length;
+    for (size_t k = 0; k < m_complex_vector_size; ++k)
+    {
+      m_gamma[k][0] /= denom;
+      m_gamma[k][1] /= denom;
+    }
+    if (length % 2 == 0)
+    {
+      m_gamma[m_complex_vector_size -1][0] /= 2;
+    }
+  }
+
+  cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old)  = delete;
+
+  cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete;
+
+  cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete;
+
+  __float128 operator()(__float128 t) const
+  {
+    using std::sin;
+    using std::cos;
+    using boost::math::constants::two_pi;
+    using std::exp;
+    __float128 s = m_gamma[0][0];
+    __float128 x = two_pi<__float128>()*(t - m_t0)/m_T;
+    fftwq_complex z;
+    z[0] = cosq(x);
+    z[1] = sinq(x);
+    fftwq_complex b{0, 0};
+    fftwq_complex u;
+    for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+    {
+      u[0] = b[0]*z[0] - b[1]*z[1];
+      u[1] = b[0]*z[1] + b[1]*z[0];
+      b[0] = m_gamma[k][0] + u[0];
+      b[1] = m_gamma[k][1] + u[1];
+    }
+
+    s += 2*(b[0]*z[0] - b[1]*z[1]);
+    return s;
+  }
+
+  __float128 prime(__float128 t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      __float128 x = two_pi<__float128>()*(t - m_t0)/m_T;
+      fftwq_complex z;
+      z[0] = cosq(x);
+      z[1] = sinq(x);
+      fftwq_complex b{0, 0};
+      // u = b*z
+      fftwq_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*m_gamma[k][0] + u[0];
+        b[1] = k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<__float128>()*(b[1]*z[0] + b[0]*z[1])/m_T;
+  }
+
+  __float128 double_prime(__float128 t) const
+  {
+      using std::sin;
+      using std::cos;
+      using boost::math::constants::two_pi;
+      using std::exp;
+      __float128 x = two_pi<__float128>()*(t - m_t0)/m_T;
+      fftwq_complex z;
+      z[0] = cosq(x);
+      z[1] = sinq(x);
+      fftwq_complex b{0, 0};
+      // u = b*z
+      fftwq_complex u;
+      for (size_t k = m_complex_vector_size - 1; k >= 1; --k)
+      {
+        u[0] = b[0]*z[0] - b[1]*z[1];
+        u[1] = b[0]*z[1] + b[1]*z[0];
+        b[0] = k*k*m_gamma[k][0] + u[0];
+        b[1] = k*k*m_gamma[k][1] + u[1];
+      }
+      // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1])
+      return -2*two_pi<__float128>()*two_pi<__float128>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T);
+  }
+
+  __float128 period() const
+  {
+    return m_T;
+  }
+
+  __float128 integrate() const
+  {
+    return m_T*m_gamma[0][0];
+  }
+
+  __float128 squared_l2() const
+  {
+    __float128 s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+      s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
+  ~cardinal_trigonometric_detail()
+  {
+    if (m_gamma)
+    {
+      fftwq_free(m_gamma);
+      m_gamma = nullptr;
+    }
+  }
+
+
+private:
+  __float128 m_t0;
+  __float128 m_h;
+  __float128 m_T;
+  fftwq_complex* m_gamma;
+  size_t m_complex_vector_size;
+};
+#endif
+
+}}}}
+#endif

boost/math/interpolators/detail/cubic_b_spline_detail.hpp

-// Copyright Nick Thompson, 2017
-// 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)
-
-#ifndef CUBIC_B_SPLINE_DETAIL_HPP
-#define CUBIC_B_SPLINE_DETAIL_HPP
-
-#include <limits>
-#include <cmath>
-#include <vector>
-#include <memory>
-#include <boost/math/constants/constants.hpp>
-#include <boost/math/special_functions/fpclassify.hpp>
-
-namespace boost{ namespace math{ namespace detail{
-
-
-template <class Real>
-class cubic_b_spline_imp
-{
-public:
-    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
-    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
-    template <class BidiIterator>
-    cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
-                       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
-                       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
-
-    Real operator()(Real x) const;
-
-    Real prime(Real x) const;
-
-private:
-    std::vector<Real> m_beta;
-    Real m_h_inv;
-    Real m_a;
-    Real m_avg;
-};
-
-
-
-template <class Real>
-Real b3_spline(Real x)
-{
-    using std::abs;
-    Real absx = abs(x);
-    if (absx < 1)
-    {
-        Real y = 2 - absx;
-        Real z = 1 - absx;
-        return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
-    }
-    if (absx < 2)
-    {
-        Real y = 2 - absx;
-        return boost::math::constants::sixth<Real>()*y*y*y;
-    }
-    return (Real) 0;
-}
-
-template<class Real>
-Real b3_spline_prime(Real x)
-{
-    if (x < 0)
-    {
-        return -b3_spline_prime(-x);
-    }
-
-    if (x < 1)
-    {
-        return x*(3*boost::math::constants::half<Real>()*x - 2);
-    }
-    if (x < 2)
-    {
-        return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
-    }
-    return (Real) 0;
-}
-
-
-template <class Real>
-template <class BidiIterator>
-cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
-                                             Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
-{
-    using boost::math::constants::third;
-
-    std::size_t length = end_p - f;
-
-    if (length < 5)
-    {
-        if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
-        {
-            throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
-        }
-        if (length < 3)
-        {
-            throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
-        }
-    }
-
-    if (boost::math::isnan(left_endpoint))
-    {
-        throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
-    }
-    if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
-    {
-        throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
-    }
-    if (step_size <= 0)
-    {
-        throw std::logic_error("The step size must be strictly > 0.\n");
-    }
-    // Storing the inverse of the stepsize does provide a measurable speedup.
-    // It's not huge, but nonetheless worthwhile.
-    m_h_inv = 1/step_size;
-
-    // Following Kress's notation, s'(a) = a1, s'(b) = b1
-    Real a1 = left_endpoint_derivative;
-    // See the finite-difference table on Wikipedia for reference on how
-    // to construct high-order estimates for one-sided derivatives:
-    // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
-    // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
-    if (boost::math::isnan(a1))
-    {
-        // For simple functions (linear, quadratic, so on)
-        // almost all the error comes from derivative estimation.
-        // This does pairwise summation which gives us another digit of accuracy over naive summation.
-        Real t0 = 4*(f[1] + third<Real>()*f[3]);
-        Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
-        a1 = m_h_inv*(t0 + t1);
-    }
-
-    Real b1 = right_endpoint_derivative;
-    if (boost::math::isnan(b1))
-    {
-        size_t n = length - 1;
-        Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
-        Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4  - 3*f[n - 2];
-
-        b1 = m_h_inv*(t0 + t1);
-    }
-
-    // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
-    // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
-    m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
-
-    // Since the splines have compact support, they decay to zero very fast outside the endpoints.
-    // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
-    // boundary [a,b] without massive error.
-    // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
-    // This algorithm for computing the average is recommended in
-    // http://www.heikohoffmann.de/htmlthesis/node134.html
-    Real t = 1;
-    for (size_t i = 0; i < length; ++i)
-    {
-        if (boost::math::isnan(f[i]))
-        {
-            std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
-            throw std::logic_error(err);
-        }
-        m_avg += (f[i] - m_avg) / t;
-        t += 1;
-    }
-
-
-    // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
-    // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
-    // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
-    // See Kress, equations 8.41
-    // The the "tridiagonal" matrix is:
-    // 1  0 -1
-    // 1  4  1
-    //    1  4  1
-    //       1  4  1
-    //          ....
-    //          1  4  1
-    //          1  0 -1
-    // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
-    std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
-    std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
-
-    rhs[0] = -2*step_size*a1;
-    rhs[rhs.size() - 1] = -2*step_size*b1;
-
-    super_diagonal[0] = 0;
-
-    for(size_t i = 1; i < rhs.size() - 1; ++i)
-    {
-        rhs[i] = 6*(f[i - 1] - m_avg);
-        super_diagonal[i] = 1;
-    }
-
-
-    // One step of row reduction on the first row to patch up the 5-diagonal problem:
-    // 1 0 -1 | r0
-    // 1 4 1  | r1
-    // mapsto:
-    // 1 0 -1 | r0
-    // 0 4 2  | r1 - r0
-    // mapsto
-    // 1 0 -1 | r0
-    // 0 1 1/2| (r1 - r0)/4
-    super_diagonal[1] = 0.5;
-    rhs[1] = (rhs[1] - rhs[0])/4;
-
-    // Now do a tridiagonal row reduction the standard way, until just before the last row:
-    for (size_t i = 2; i < rhs.size() - 1; ++i)
-    {
-        Real diagonal = 4 - super_diagonal[i - 1];
-        rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
-        super_diagonal[i] /= diagonal;
-    }
-
-    // Now the last row, which is in the form
-    // 1 sd[n-3] 0      | rhs[n-3]
-    // 0  1     sd[n-2] | rhs[n-2]
-    // 1  0     -1      | rhs[n-1]
-    Real final_subdiag = -super_diagonal[rhs.size() - 3];
-    rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
-    Real final_diag = -1/final_subdiag;
-    // Now we're here:
-    // 1 sd[n-3] 0         | rhs[n-3]
-    // 0  1     sd[n-2]    | rhs[n-2]
-    // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
-
-    final_diag = final_diag - super_diagonal[rhs.size() - 2];
-    rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
-
-
-    // Back substitutions:
-    m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
-    for(size_t i = rhs.size() - 2; i > 0; --i)
-    {
-        m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
-    }
-    m_beta[0] = m_beta[2] + rhs[0];
-}
-
-template<class Real>
-Real cubic_b_spline_imp<Real>::operator()(Real x) const
-{
-    // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
-    // just the (at most 5) whose support overlaps the argument.
-    Real z = m_avg;
-    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) (max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0);
-    for (size_t k = k_min; k <= k_max; ++k)
-    {
-        z += m_beta[k]*b3_spline(t - k);
-    }
-
-    return z;
-}
-
-template<class Real>
-Real cubic_b_spline_imp<Real>::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_prime(t - k);
-    }
-    return z*m_h_inv;
-}
-
-}}}
-#endif
+// Copyright Nick Thompson, 2017
+// 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)
+
+#ifndef CUBIC_B_SPLINE_DETAIL_HPP
+#define CUBIC_B_SPLINE_DETAIL_HPP
+
+#include <limits>
+#include <cmath>
+#include <vector>
+#include <memory>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+
+template <class Real>
+class cubic_b_spline_imp
+{
+public:
+    // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
+    // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
+    template <class BidiIterator>
+    cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                       Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+                       Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
+
+    Real operator()(Real x) const;
+
+    Real prime(Real x) const;
+
+    Real double_prime(Real x) const;
+
+private:
+    std::vector<Real> m_beta;
+    Real m_h_inv;
+    Real m_a;
+    Real m_avg;
+};
+
+
+
+template <class Real>
+Real b3_spline(Real x)
+{
+    using std::abs;
+    Real absx = abs(x);
+    if (absx < 1)
+    {
+        Real y = 2 - absx;
+        Real z = 1 - absx;
+        return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
+    }
+    if (absx < 2)
+    {
+        Real y = 2 - absx;
+        return boost::math::constants::sixth<Real>()*y*y*y;
+    }
+    return (Real) 0;
+}
+
+template<class Real>
+Real b3_spline_prime(Real x)
+{
+    if (x < 0)
+    {
+        return -b3_spline_prime(-x);
+    }
+
+    if (x < 1)
+    {
+        return x*(3*boost::math::constants::half<Real>()*x - 2);
+    }
+    if (x < 2)
+    {
+        return -boost::math::constants::half<Real>()*(2 - x)*(2 - 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*x - 2;
+    }
+    if (x < 2)
+    {
+        return (2 - x);
+    }
+    return (Real) 0;
+}
+
+
+template <class Real>
+template <class BidiIterator>
+cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
+                                             Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
+{
+    using boost::math::constants::third;
+
+    std::size_t length = end_p - f;
+
+    if (length < 5)
+    {
+        if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
+        {
+            throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
+        }
+        if (length < 3)
+        {
+            throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
+        }
+    }
+
+    if (boost::math::isnan(left_endpoint))
+    {
+        throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
+    }
+    if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
+    {
+        throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
+    }
+    if (step_size <= 0)
+    {
+        throw std::logic_error("The step size must be strictly > 0.\n");
+    }
+    // Storing the inverse of the stepsize does provide a measurable speedup.
+    // It's not huge, but nonetheless worthwhile.
+    m_h_inv = 1/step_size;
+
+    // Following Kress's notation, s'(a) = a1, s'(b) = b1
+    Real a1 = left_endpoint_derivative;
+    // See the finite-difference table on Wikipedia for reference on how
+    // to construct high-order estimates for one-sided derivatives:
+    // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
+    // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
+    if (boost::math::isnan(a1))
+    {
+        // For simple functions (linear, quadratic, so on)
+        // almost all the error comes from derivative estimation.
+        // This does pairwise summation which gives us another digit of accuracy over naive summation.
+        Real t0 = 4*(f[1] + third<Real>()*f[3]);
+        Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
+        a1 = m_h_inv*(t0 + t1);
+    }
+
+    Real b1 = right_endpoint_derivative;
+    if (boost::math::isnan(b1))
+    {
+        size_t n = length - 1;
+        Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
+        Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4  - 3*f[n - 2];
+
+        b1 = m_h_inv*(t0 + t1);
+    }
+
+    // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
+    // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
+    m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+    // Since the splines have compact support, they decay to zero very fast outside the endpoints.
+    // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
+    // boundary [a,b] without massive error.
+    // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
+    // This algorithm for computing the average is recommended in
+    // http://www.heikohoffmann.de/htmlthesis/node134.html
+    Real t = 1;
+    for (size_t i = 0; i < length; ++i)
+    {
+        if (boost::math::isnan(f[i]))
+        {
+            std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
+            throw std::logic_error(err);
+        }
+        m_avg += (f[i] - m_avg) / t;
+        t += 1;
+    }
+
+
+    // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
+    // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
+    // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
+    // See Kress, equations 8.41
+    // The the "tridiagonal" matrix is:
+    // 1  0 -1
+    // 1  4  1
+    //    1  4  1
+    //       1  4  1
+    //          ....
+    //          1  4  1
+    //          1  0 -1
+    // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
+    std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
+    std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
+
+    rhs[0] = -2*step_size*a1;
+    rhs[rhs.size() - 1] = -2*step_size*b1;
+
+    super_diagonal[0] = 0;
+
+    for(size_t i = 1; i < rhs.size() - 1; ++i)
+    {
+        rhs[i] = 6*(f[i - 1] - m_avg);
+        super_diagonal[i] = 1;
+    }
+
+
+    // One step of row reduction on the first row to patch up the 5-diagonal problem:
+    // 1 0 -1 | r0
+    // 1 4 1  | r1
+    // mapsto:
+    // 1 0 -1 | r0
+    // 0 4 2  | r1 - r0
+    // mapsto
+    // 1 0 -1 | r0
+    // 0 1 1/2| (r1 - r0)/4
+    super_diagonal[1] = 0.5;
+    rhs[1] = (rhs[1] - rhs[0])/4;
+
+    // Now do a tridiagonal row reduction the standard way, until just before the last row:
+    for (size_t i = 2; i < rhs.size() - 1; ++i)
+    {
+        Real diagonal = 4 - super_diagonal[i - 1];
+        rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
+        super_diagonal[i] /= diagonal;
+    }
+
+    // Now the last row, which is in the form
+    // 1 sd[n-3] 0      | rhs[n-3]
+    // 0  1     sd[n-2] | rhs[n-2]
+    // 1  0     -1      | rhs[n-1]
+    Real final_subdiag = -super_diagonal[rhs.size() - 3];
+    rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
+    Real final_diag = -1/final_subdiag;
+    // Now we're here:
+    // 1 sd[n-3] 0         | rhs[n-3]
+    // 0  1     sd[n-2]    | rhs[n-2]
+    // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
+
+    final_diag = final_diag - super_diagonal[rhs.size() - 2];
+    rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
+
+
+    // Back substitutions:
+    m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
+    for(size_t i = rhs.size() - 2; i > 0; --i)
+    {
+        m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
+    }
+    m_beta[0] = m_beta[2] + rhs[0];
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::operator()(Real x) const
+{
+    // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
+    // just the (at most 5) whose support overlaps the argument.
+    Real z = m_avg;
+    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) (max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0);
+    for (size_t k = k_min; k <= k_max; ++k)
+    {
+        z += m_beta[k]*b3_spline(t - k);
+    }
+
+    return z;
+}
+
+template<class Real>
+Real cubic_b_spline_imp<Real>::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_prime(t - k);
+    }
+    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

boost/math/interpolators/detail/cubic_hermite_detail.hpp

+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
+#include <stdexcept>
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <sstream>
+#include <limits>
+
+namespace boost::math::interpolators::detail {
+
+template<class RandomAccessContainer>
+class cubic_hermite_detail {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    cubic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer dydx)
+     : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}
+    {
+        using std::abs;
+        using std::isnan;
+        if (x_.size() != y_.size())
+        {
+            throw std::domain_error("There must be the same number of ordinates as abscissas.");
+        }
+        if (x_.size() != dydx_.size())
+        {
+            throw std::domain_error("There must be the same number of ordinates as derivative values.");
+        }
+        if (x_.size() < 2)
+        {
+            throw std::domain_error("Must be at least two data points.");
+        }
+        Real x0 = x_[0];
+        for (size_t i = 1; i < x_.size(); ++i)
+        {
+            Real x1 = x_[i];
+            if (x1 <= x0)
+            {
+                std::ostringstream oss;
+                oss.precision(std::numeric_limits<Real>::digits10+3);
+                oss << "Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}, ";
+                oss << "but at x[" << i - 1 << "] = " << x0 << ", and x[" << i << "] = " << x1 << ".\n";
+                throw std::domain_error(oss.str());
+            }
+            x0 = x1;
+        }
+    }
+
+    void push_back(Real x, Real y, Real dydx)
+    {
+        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);
+    }
+
+    Real operator()(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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 s0 = dydx_[i];
+        Real s1 = dydx_[i+1];
+        Real dx = (x1-x0);
+        Real t = (x-x0)/dx;
+
+        // See the section 'Representations' in the page
+        // https://en.wikipedia.org/wiki/Cubic_Hermite_spline
+        Real y = (1-t)*(1-t)*(y0*(1+2*t) + s0*(x-x0))
+              + t*t*(y1*(3-2*t) + dx*s1*(t-1));
+        return y;
+    }
+
+    Real prime(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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 y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real s0 = dydx_[i];
+        Real s1 = dydx_[i+1];
+        Real dx = (x1-x0);
+
+        Real d1 = (y1 - y0 - s0*dx)/(dx*dx);
+        Real d2 = (s1 - s0)/(2*dx);
+        Real c2 = 3*d1 - 2*d2;
+        Real c3 = 2*(d2 - d1)/dx;
+        return s0 + 2*c2*(x-x0) + 3*c3*(x-x0)*(x-x0); 
+    }
+
+
+    friend std::ostream& operator<<(std::ostream & os, const cubic_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] << "),  ";
+        }
+        auto n = m.x_.size()-1;
+        os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ")}";
+        return os;
+    }
+
+    auto size() const
+    {
+        return x_.size();
+    }
+
+    int64_t bytes() const
+    {
+        return 3*x_.size()*sizeof(Real) + 3*sizeof(x_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return {x_.front(), x_.back()};
+    }
+
+    RandomAccessContainer x_;
+    RandomAccessContainer y_;
+    RandomAccessContainer dydx_;
+};
+
+template<class RandomAccessContainer>
+class cardinal_cubic_hermite_detail {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    cardinal_cubic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer dydx, Real x0, Real dx)
+    : y_{std::move(y)}, dy_{std::move(dydx)}, x0_{x0}, inv_dx_{1/dx}
+    {
+        using std::abs;
+        using std::isnan;
+        if (y_.size() != dy_.size())
+        {
+            throw std::domain_error("There must be the same number of derivatives as ordinates.");
+        }
+        if (y_.size() < 2)
+        {
+            throw std::domain_error("Must be at least two data points.");
+        }
+        if (dx <= 0)
+        {
+            throw std::domain_error("dx > 0 is required.");
+        }
+
+        for (auto & dy : dy_)
+        {
+            dy *= dx;
+        }
+    }
+
+    // Why not implement push_back? It's awkward: If the buffer is circular, x0_ += dx_.
+    // If the buffer is not circular, x0_ is unchanged.
+    // We need a concept for circular_buffer!
+
+    inline Real operator()(Real x) const
+    {
+        const Real xf = x0_ + (y_.size()-1)/inv_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 y_.back();
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(y_.size())>(ii);
+        Real t = s - ii;
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+
+        Real r = 1-t;
+        return r*r*(y0*(1+2*t) + dy0*t)
+              + t*t*(y1*(3-2*t) - dy1*r);
+    }
+
+    inline Real prime(Real x) const
+    {
+        const Real xf = x0_ + (y_.size()-1)/inv_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 dy_.back()*inv_dx_;
+        }
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(y_.size())>(ii);
+        Real t = s - ii;
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+
+        Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
+        return dy*inv_dx_;
+    }
+
+
+    auto size() const
+    {
+        return y_.size();
+    }
+
+    int64_t bytes() const
+    {
+        return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
+private:
+
+    RandomAccessContainer y_;
+    RandomAccessContainer dy_;
+    Real x0_;
+    Real inv_dx_;
+};
+
+
+template<class RandomAccessContainer>
+class cardinal_cubic_hermite_detail_aos {
+public:
+    using Point = typename RandomAccessContainer::value_type;
+    using Real = typename Point::value_type;
+
+    cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
+    : dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
+    {
+        if (dat_.size() < 2)
+        {
+            throw std::domain_error("Must be at least two data points.");
+        }
+        if (dat_[0].size() != 2)
+        {
+            throw std::domain_error("Each datum must contain (y, y'), and nothing else.");
+        }
+        if (dx <= 0)
+        {
+            throw std::domain_error("dx > 0 is required.");
+        }
+
+        for (auto & d : dat_)
+        {
+            d[1] *= dx;
+        }
+    }
+
+    inline Real operator()(Real x) const
+    {
+        const Real xf = x0_ + (dat_.size()-1)/inv_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 dat_.back()[0];
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(dat_.size())>(ii);
+
+        Real t = s - ii;
+        // If we had infinite precision, this would never happen.
+        // But we don't have infinite precision.
+        if (t == 0)
+        {
+            return dat_[i][0];
+        }
+        Real y0 = dat_[i][0];
+        Real y1 = dat_[i+1][0];
+        Real dy0 = dat_[i][1];
+        Real dy1 = dat_[i+1][1];
+
+        Real r = 1-t;
+        return r*r*(y0*(1+2*t) + dy0*t)
+              + t*t*(y1*(3-2*t) - dy1*r);
+    }
+
+    inline Real prime(Real x) const
+    {
+        const Real xf = x0_ + (dat_.size()-1)/inv_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 dat_.back()[1]*inv_dx_;
+        }
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(dat_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0)
+        {
+            return dat_[i][1]*inv_dx_;
+        }
+        Real y0 = dat_[i][0];
+        Real dy0 = dat_[i][1];
+        Real y1 = dat_[i+1][0];
+        Real dy1 = dat_[i+1][1];
+
+        Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
+        return dy*inv_dx_;
+    }
+
+
+    auto size() const
+    {
+        return dat_.size();
+    }
+
+    int64_t bytes() const
+    {
+        return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (dat_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
+
+private:
+    RandomAccessContainer dat_;
+    Real x0_;
+    Real inv_dx_;
+};
+
+
+}
+#endif

boost/math/interpolators/detail/quintic_hermite_detail.hpp

+/*
+ * 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)
+ */
+#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
+#include <algorithm>
+#include <stdexcept>
+#include <sstream>
+#include <cmath>
+
+namespace boost::math::interpolators::detail {
+
+template<class RandomAccessContainer>
+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);
+    }
+
+    inline Real operator()(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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;
+        Real t2 = t*t;
+        Real t3 = t2*t;
+
+        // 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- t3*(10 + t*(-15 + 6*t)))*y0;
+        y += t*(1+ t2*(-6 + t*(8 -3*t)))*v0*dx;
+        y += t2*(1 + t*(-3 + t*(3-t)))*a0*dx*dx/2;
+        y += t3*((1 + t*(-2 + t))*a1*dx*dx/2 + (-4 + t*(7 - 3*t))*v1*dx + (10 + t*(-15 + 6*t))*y1);
+        return y;
+    }
+
+    inline Real prime(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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 dx = x1 - x0;
+
+        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 t= (x-x0)/dx;
+        Real t2 = t*t;
+
+        Real dydx = 30*t2*(1 - 2*t + t*t)*(y1-y0)/dx;
+        dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
+        dydx += (t*dx/2)*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
+        return dydx;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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 d2ydx2_.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 dx = x1 - x0;
+
+        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 t = (x-x0)/dx;
+
+        Real d2ydx2 = 60*t*(1 + t*(-3 + 2*t))*(y1-y0)/(dx*dx);
+        d2ydx2 += 12*t*(-3 + t*(8 - 5*t))*v0/dx;
+        d2ydx2 -= 12*t*(2 + t*(-7 + 5*t))*v1/dx;
+        d2ydx2 += (1 + t*(-9 + t*(18 - 10*t)))*a0;
+        d2ydx2 += t*(3 + t*(-12 + 10*t))*a1;
+
+        return d2ydx2;
+    }
+
+    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;
+    }
+
+    int64_t bytes() const
+    {
+        return 4*x_.size()*sizeof(x_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return {x_.front(), x_.back()};
+    }
+
+private:
+    RandomAccessContainer x_;
+    RandomAccessContainer y_;
+    RandomAccessContainer dydx_;
+    RandomAccessContainer d2ydx2_;
+};
+
+
+template<class RandomAccessContainer>
+class cardinal_quintic_hermite_detail {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    cardinal_quintic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, Real x0, Real dx)
+    : y_{std::move(y)}, dy_{std::move(dydx)}, d2y_{std::move(d2ydx2)}, x0_{x0}, inv_dx_{1/dx}
+    {
+        if (y_.size() != dy_.size())
+        {
+            throw std::domain_error("Numbers of derivatives must = number of abscissas.");
+        }
+        if (y_.size() != d2y_.size())
+        {
+            throw std::domain_error("Number of second derivatives must equal number of abscissas.");
+        }
+        if (y_.size() < 2)
+        {
+            throw std::domain_error("At least 2 abscissas are required.");
+        }
+        if (dx <= 0)
+        {
+            throw std::domain_error("dx > 0 is required.");
+        }
+
+        for (auto & dy : dy_)
+        {
+            dy *= dx;
+        }
+
+        for (auto & d2y : d2y_)
+        {
+            d2y *= (dx*dx)/2;
+        }
+    }
+
+
+    inline Real operator()(Real x) const
+    {
+        const Real xf = x0_ + (y_.size()-1)/inv_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 y_.back();
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(y_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0)
+        {
+            return y_[i];
+        }
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real d2y0 = d2y_[i];
+        Real d2y1 = d2y_[i+1];
+
+        // 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 + t*(-15 + 6*t)))*y0;
+        y += t*(1+ t*t*(-6 + t*(8 -3*t)))*dy0;
+        y += t*t*(1 + t*(-3 + t*(3-t)))*d2y0;
+        y += t*t*t*((1 + t*(-2 + t))*d2y1 + (-4 + t*(7 -3*t))*dy1 + (10 + t*(-15 + 6*t))*y1);
+        return y;
+    }
+
+    inline Real prime(Real x) const
+    {
+        const Real xf = x0_ + (y_.size()-1)/inv_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 dy_.back()*inv_dx_;
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(y_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0)
+        {
+            return dy_[i]*inv_dx_;
+        }
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real d2y0 = d2y_[i];
+        Real d2y1 = d2y_[i+1];
+
+        Real dydx = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
+        dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*dy0 - t*t*(12 - 28*t + 15*t*t)*dy1;
+        dydx += t*((2 - 9*t + 12*t*t - 5*t*t*t)*d2y0 + t*(3 - 8*t + 5*t*t)*d2y1);
+        dydx *= inv_dx_;
+        return dydx;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        const Real xf = x0_ + (y_.size()-1)/inv_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 d2y_.back()*2*inv_dx_*inv_dx_;
+        }
+
+        return this->unchecked_double_prime(x);
+    }
+
+    inline Real unchecked_double_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(y_.size())>(ii);
+        Real t = s - ii;
+        if (t==0)
+        {
+            return d2y_[i]*2*inv_dx_*inv_dx_;
+        }
+
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real d2y0 = d2y_[i];
+        Real d2y1 = d2y_[i+1];
+
+        Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
+        d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
+        d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
+        return d2ydx2;
+    }
+
+    int64_t bytes() const
+    {
+        return 3*y_.size()*sizeof(Real) + 2*sizeof(Real);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
+private:
+    RandomAccessContainer y_;
+    RandomAccessContainer dy_;
+    RandomAccessContainer d2y_;
+    Real x0_;
+    Real inv_dx_;
+};
+
+
+template<class RandomAccessContainer>
+class cardinal_quintic_hermite_detail_aos {
+public:
+    using Point = typename RandomAccessContainer::value_type;
+    using Real = typename Point::value_type;
+    cardinal_quintic_hermite_detail_aos(RandomAccessContainer && data, Real x0, Real dx)
+    : data_{std::move(data)} , x0_{x0}, inv_dx_{1/dx}
+    {
+        if (data_.size() < 2)
+        {
+            throw std::domain_error("At least two points are required to interpolate using cardinal_quintic_hermite_aos");
+        }
+
+        if (data_[0].size() != 3)
+        {
+            throw std::domain_error("Each datum passed to the cardinal_quintic_hermite_aos must have three elements: {y, y', y''}");
+        }
+        if (dx <= 0)
+        {
+            throw std::domain_error("dx > 0 is required.");
+        }
+
+        for (auto & datum : data_)
+        {
+            datum[1] *= dx;
+            datum[2] *= (dx*dx/2);
+        }
+    }
+
+
+    inline Real operator()(Real x) const
+    {
+        const Real xf = x0_ + (data_.size()-1)/inv_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()[0];
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(data_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0)
+        {
+            return data_[i][0];
+        }
+
+        Real y0 = data_[i][0];
+        Real dy0 = data_[i][1];
+        Real d2y0 = data_[i][2];
+        Real y1 = data_[i+1][0];
+        Real dy1 = data_[i+1][1];
+        Real d2y1 = data_[i+1][2];
+
+        Real y = (1 - t*t*t*(10 + t*(-15 + 6*t)))*y0;
+        y += t*(1 + t*t*(-6 + t*(8 - 3*t)))*dy0;
+        y += t*t*(1 + t*(-3 + t*(3 - t)))*d2y0;
+        y += t*t*t*((1 + t*(-2 + t))*d2y1 + (-4 + t*(7 - 3*t))*dy1 + (10 + t*(-15 + 6*t))*y1);
+        return y;
+    }
+
+    inline Real prime(Real x) const
+    {
+        const Real xf = x0_ + (data_.size()-1)/inv_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]*inv_dx_;
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(data_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0)
+        {
+            return data_[i][1]*inv_dx_;
+        }
+
+
+        Real y0 = data_[i][0];
+        Real y1 = data_[i+1][0];
+        Real v0 = data_[i][1];
+        Real v1 = data_[i+1][1];
+        Real a0 = data_[i][2];
+        Real a1 = data_[i+1][2];
+
+        Real dy = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
+        dy += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
+        dy += t*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
+        return dy*inv_dx_;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        const Real xf = x0_ + (data_.size()-1)/inv_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()[2]*2*inv_dx_*inv_dx_;
+        }
+
+        return this->unchecked_double_prime(x);
+    }
+
+    inline Real unchecked_double_prime(Real x) const
+    {
+        using std::floor;
+        Real s = (x-x0_)*inv_dx_;
+        Real ii = floor(s);
+        auto i = static_cast<decltype(data_.size())>(ii);
+        Real t = s - ii;
+        if (t == 0) {
+            return data_[i][2]*2*inv_dx_*inv_dx_;
+        }
+        Real y0 = data_[i][0];
+        Real dy0 = data_[i][1];
+        Real d2y0 = data_[i][2];
+        Real y1 = data_[i+1][0];
+        Real dy1 = data_[i+1][1];
+        Real d2y1 = data_[i+1][2];
+
+        Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
+        d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
+        d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
+        return d2ydx2;
+    }
+
+    int64_t bytes() const
+    {
+        return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (data_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
+private:
+    RandomAccessContainer data_;
+    Real x0_;
+    Real inv_dx_;
+};
+
+}
+#endif

boost/math/interpolators/detail/septic_hermite_detail.hpp

+/*
+ * 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)
+ */
+#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP
+#include <algorithm>
+#include <stdexcept>
+#include <sstream>
+#include <cmath>
+
+namespace boost::math::interpolators::detail {
+
+template<class RandomAccessContainer>
+class septic_hermite_detail {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    septic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3) 
+    : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)}
+    {
+        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() != d3ydx3_.size())
+        {
+            throw std::domain_error("Number of third 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, Real d3ydx3)
+    {
+        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);
+        d3ydx3_.push_back(d3ydx3);
+    }
+
+    Real operator()(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::digits10+3);
+            oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
+                << x_[0] << ", " << x_.back() << "]";
+            throw std::domain_error(oss.str());
+        }
+        // t \in [0, 1)
+        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 dx = (x1-x0);
+        Real t = (x-x0)/dx;
+
+        // 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 dx2 = dx*dx/2;
+        Real dx3 = dx2*dx/3;
+
+        Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
+        Real z4 = -s;
+        Real z0 = s + 1;
+        Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36 + 10*t))));
+        Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15 + 4*t))));
+        Real z3 = t3*(1 + t*(-4 + t*(6 + t*(-4 + t))));
+        Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
+        Real z6 = t4*(5 + t*(-14 + t*(13 - 4*t)));
+        Real z7 = t4*(-1 + t*(3 + t*(-3+t)));
+
+        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];
+
+        return z0*y0 + z4*y1 + (z1*v0 + z5*v1)*dx + (z2*a0 + z6*a1)*dx2 + (z3*j0 + z7*j1)*dx3;
+    }
+
+    Real prime(Real x) const
+    {
+        if  (x < x_[0] || x > x_.back())
+        {
+            std::ostringstream oss;
+            oss.precision(std::numeric_limits<Real>::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 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 j0 = d3ydx3_[i];
+        Real j1 = d3ydx3_[i+1];
+        Real dx = x1 - x0;
+        Real t = (x-x0)/dx;
+        Real t2 = t*t;
+        Real t3 = t2*t;
+        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
+        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
+        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
+        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
+        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
+        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
+        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
+
+        Real dydx = z0*(y1-y0)/dx;
+        dydx += z1*v0 + z2*v1;
+        dydx += (x-x0)*(z3*a0 + z4*a1);
+        dydx += (x-x0)*(x-x0)*(z5*j0 + z6*j1)/6;
+        return dydx;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+
+    friend std::ostream& operator<<(std::ostream & os, const septic_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] <<  ", " << m.d3ydx3_[i] << "),  ";
+        }
+        auto n = m.x_.size()-1;
+        os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << m.d3ydx3_[n] << ")}";
+        return os;
+    }
+
+    int64_t bytes()
+    {
+        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_;
+    RandomAccessContainer dydx_;
+    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)}, dy_{std::move(dydx)}, d2y_{std::move(d2ydx2)}, d3y_{std::move(d3ydx3)}, x0_{x0}, inv_dx_{1/dx}
+    {
+        if (y_.size() != dy_.size())
+        {
+            throw std::domain_error("Numbers of derivatives must = number of ordinates.");
+        }
+        if (y_.size() != d2y_.size())
+        {
+            throw std::domain_error("Number of second derivatives must equal number of ordinates.");
+        }
+        if (y_.size() != d3y_.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.");
+        }
+
+        if (dx <= 0)
+        {
+            throw std::domain_error("dx > 0 is required.");
+        }
+
+        for (auto & dy : dy_)
+        {
+            dy *= dx;
+        }
+        for (auto & d2y : d2y_)
+        {
+            d2y *= (dx*dx/2);
+        }
+        for (auto & d3y : d3y_)
+        {
+            d3y *= (dx*dx*dx/6);
+        }
+
+    }
+
+    inline Real operator()(Real x) const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_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 y_.back();
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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;
+        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*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
+        Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
+        Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
+        Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
+        Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
+        Real z7 = t4*(-1 + t*(3+t*(-3+t)));
+
+        Real y0 = y_[i];
+        Real y1 = y_[i+1];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real a0 = d2y_[i];
+        Real a1 = d2y_[i+1];
+        Real j0 = d3y_[i];
+        Real j1 = d3y_[i+1];
+
+        return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1;
+    }
+
+    inline Real prime(Real x) const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_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 dy_.back()/inv_dx_;
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real a0 = d2y_[i];
+        Real a1 = d2y_[i+1];
+        Real j0 = d3y_[i];
+        Real j1 = d3y_[i+1];
+        Real t2 = t*t;
+        Real t3 = t2*t;
+        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
+        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
+        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
+        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
+        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
+        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
+        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
+
+        Real dydx = z0*(y1-y0)*inv_dx_;
+        dydx += (z1*dy0 + z2*dy1)*inv_dx_;
+        dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
+        dydx += t*t*(z5*j0 + z6*j1);
+        return dydx;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_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 d2y_.back()*2*inv_dx_*inv_dx_;
+        }
+
+        return this->unchecked_double_prime(x);
+    }
+
+    inline Real unchecked_double_prime(Real x) const
+    {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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];
+        Real dy0 = dy_[i];
+        Real dy1 = dy_[i+1];
+        Real a0 = d2y_[i];
+        Real a1 = d2y_[i+1];
+        Real j0 = d3y_[i];
+        Real j1 = d3y_[i+1];
+        Real t2 = t*t;
+
+        Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
+        Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
+        Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
+        Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
+        Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
+        Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
+        Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
+
+        Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
+        d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
+        d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
+        d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
+
+        return d2ydx2;
+    }
+
+    int64_t bytes() const
+    {
+        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_;
+    RandomAccessContainer d2y_;
+    RandomAccessContainer d3y_;
+    Real x0_;
+    Real inv_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}, inv_dx_{1/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.");
+        }
+
+        for (auto & datum : data_)
+        {
+            datum[1] *= dx;
+            datum[2] *= (dx*dx/2);
+            datum[3] *= (dx*dx*dx/6);
+        }
+    }
+
+    inline Real operator()(Real x) const
+    {
+        Real xf = x0_ + (data_.size()-1)/inv_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()[0];
+        }
+        return this->unchecked_evaluation(x);
+    }
+
+    inline Real unchecked_evaluation(Real x) const
+    {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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;
+
+        Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
+        Real z4 = -s;
+        Real z0 = s + 1;
+        Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
+        Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
+        Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
+        Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
+        Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
+        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;
+    }
+
+    inline Real prime(Real x) const
+    {
+        Real xf = x0_ + (data_.size()-1)/inv_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]*inv_dx_;
+        }
+
+        return this->unchecked_prime(x);
+    }
+
+    inline Real unchecked_prime(Real x) const
+    {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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];
+        Real dy0 = data_[i][1];
+        Real dy1 = data_[i+1][1];
+        Real a0 = data_[i][2];
+        Real a1 = data_[i+1][2];
+        Real j0 = data_[i][3];
+        Real j1 = data_[i+1][3];
+        Real t2 = t*t;
+        Real t3 = t2*t;
+        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
+        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
+        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
+        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
+        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
+        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
+        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
+
+        Real dydx = z0*(y1-y0)*inv_dx_;
+        dydx += (z1*dy0 + z2*dy1)*inv_dx_;
+        dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
+        dydx += t*t*(z5*j0 + z6*j1);
+        return dydx;
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        Real xf = x0_ + (data_.size()-1)/inv_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()[2]*2*inv_dx_*inv_dx_;
+        }
+
+        return this->unchecked_double_prime(x);
+    }
+
+    inline Real unchecked_double_prime(Real x) const
+    {
+        using std::floor;
+        Real s3 = (x-x0_)*inv_dx_;
+        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];
+        Real dy1 = data_[i+1][1];
+        Real a0 = data_[i][2];
+        Real a1 = data_[i+1][2];
+        Real j0 = data_[i][3];
+        Real j1 = data_[i+1][3];
+        Real t2 = t*t;
+
+        Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
+        Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
+        Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
+        Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
+        Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
+        Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
+        Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
+
+        Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
+        d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
+        d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
+        d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
+
+        return d2ydx2;
+    }
+
+    int64_t bytes() const
+    {
+        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_;
+    Real inv_dx_;
+};
+
+
+}
+#endif

boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp

+/*
+ *  Copyright Nick Thompson, 2019
+ *  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)
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
+
+#include <vector>
+#include <utility> // for std::move
+#include <limits>
+#include <boost/assert.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class TimeContainer, class SpaceContainer>
+class vector_barycentric_rational_imp
+{
+public:
+    using Real = typename TimeContainer::value_type;
+    using Point = typename SpaceContainer::value_type;
+
+    vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order);
+
+    void operator()(Point& p, Real t) const;
+
+    void eval_with_prime(Point& x, Point& dxdt, Real t) const;
+
+    // The barycentric weights are only interesting to the unit tests:
+    Real weight(size_t i) const { return w_[i]; }
+
+private:
+
+    void calculate_weights(size_t approximation_order);
+
+    TimeContainer t_;
+    SpaceContainer y_;
+    TimeContainer w_;
+};
+
+template <class TimeContainer, class SpaceContainer>
+vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order)
+{
+    using std::numeric_limits;
+    t_ = std::move(t);
+    y_ = std::move(y);
+
+    BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points.");
+    BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length.");
+    for (size_t i = 1; i < t_.size(); ++i)
+    {
+        BOOST_ASSERT_MSG(t_[i] - t_[i-1] >  (numeric_limits<typename TimeContainer::value_type>::min)(), "The abscissas must be listed in strictly increasing order t[0] < t[1] < ... < t[n-1].");
+    }
+    calculate_weights(approximation_order);
+}
+
+
+template<class TimeContainer, class SpaceContainer>
+void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::calculate_weights(size_t approximation_order)
+{
+    using Real = typename TimeContainer::value_type;
+    using std::abs;
+    int64_t n = t_.size();
+    w_.resize(n, Real(0));
+    for(int64_t k = 0; k < n; ++k)
+    {
+        int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
+        int64_t i_max = k;
+        if (k >= n - (std::ptrdiff_t)approximation_order)
+        {
+            i_max = n - approximation_order - 1;
+        }
+
+        for(int64_t i = i_min; i <= i_max; ++i)
+        {
+            Real inv_product = 1;
+            int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
+            for(int64_t j = i; j <= j_max; ++j)
+            {
+                if (j == k)
+                {
+                    continue;
+                }
+                Real diff = t_[k] - t_[j];
+                inv_product *= diff;
+            }
+            if (i % 2 == 0)
+            {
+                w_[k] += 1/inv_product;
+            }
+            else
+            {
+                w_[k] -= 1/inv_product;
+            }
+        }
+    }
+}
+
+
+template<class TimeContainer, class SpaceContainer>
+void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const
+{
+    using Real = typename TimeContainer::value_type;
+    for (auto & x : p)
+    {
+        x = Real(0);
+    }
+    Real denominator = 0;
+    for(size_t i = 0; i < t_.size(); ++i)
+    {
+        // See associated commentary in the scalar version of this function.
+        if (t == t_[i])
+        {
+            p = y_[i];
+            return;
+        }
+        Real x = w_[i]/(t - t_[i]);
+        for (decltype(p.size()) j = 0; j < p.size(); ++j)
+        {
+            p[j] += x*y_[i][j];
+        }
+        denominator += x;
+    }
+    for (decltype(p.size()) j = 0; j < p.size(); ++j)
+    {
+        p[j] /= denominator;
+    }
+    return;
+}
+
+template<class TimeContainer, class SpaceContainer>
+void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::eval_with_prime(typename SpaceContainer::value_type& x, typename SpaceContainer::value_type& dxdt, typename TimeContainer::value_type t) const
+{
+    using Point = typename SpaceContainer::value_type;
+    using Real = typename TimeContainer::value_type;
+    this->operator()(x, t);
+    Point numerator;
+    for (decltype(x.size()) i = 0; i < x.size(); ++i)
+    {
+        numerator[i] = 0;
+    }
+    Real denominator = 0;
+    for(decltype(t_.size()) i = 0; i < t_.size(); ++i)
+    {
+        if (t == t_[i])
+        {
+            Point sum;
+            for (decltype(x.size()) i = 0; i < x.size(); ++i)
+            {
+                sum[i] = 0;
+            }
+
+            for (decltype(t_.size()) j = 0; j < t_.size(); ++j)
+            {
+                if (j == i)
+                {
+                    continue;
+                }
+                for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
+                {
+                    sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]);
+                }
+            }
+            for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
+            {
+                dxdt[k] = -sum[k]/w_[i];
+            }
+            return;
+        }
+        Real tw = w_[i]/(t - t_[i]);
+        Point diff;
+        for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
+        {
+            diff[j] = (x[j] - y_[i][j])/(t-t_[i]);
+        }
+        for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
+        {
+            numerator[j] += tw*diff[j];
+        }
+        denominator += tw;
+    }
+
+    for (decltype(dxdt.size()) j = 0; j < dxdt.size(); ++j)
+    {
+        dxdt[j] = numerator[j]/denominator;
+    }
+    return;
+}
+
+}}}
+#endif

boost/math/interpolators/detail/whittaker_shannon_detail.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP
+#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP
+#include <boost/assert.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/sin_pi.hpp>
+#include <boost/math/special_functions/cos_pi.hpp>
+
+namespace boost { namespace math { namespace interpolators { namespace detail {
+
+template<class RandomAccessContainer>
+class whittaker_shannon_detail {
+public:
+
+    using Real = typename RandomAccessContainer::value_type;
+    whittaker_shannon_detail(RandomAccessContainer&& y, Real const & t0, Real const & h) : m_y{std::move(y)}, m_t0{t0}, m_h{h}
+    {
+        for (size_t i = 1; i < m_y.size(); i += 2)
+        {
+            m_y[i] = -m_y[i];
+        }
+    }
+
+    inline Real operator()(Real t) const {
+        using boost::math::constants::pi;
+        using std::isfinite;
+        using std::floor;
+        Real y = 0;
+        Real x = (t - m_t0)/m_h;
+        Real z = x;
+        auto it = m_y.begin();
+
+        // For some reason, neither clang nor g++ will cache the address of m_y.end() in a register.
+        // Hence make a copy of it:
+        auto end = m_y.end();
+        while(it != end)
+        {
+
+            y += *it++/z;
+            z -= 1;
+        }
+
+        if (!isfinite(y))
+        {
+            BOOST_ASSERT_MSG(floor(x) == ceil(x), "Floor and ceiling should be equal.\n");
+            size_t i = static_cast<size_t>(floor(x));
+            if (i & 1)
+            {
+                return -m_y[i];
+            }
+            return m_y[i];
+        }
+        return y*boost::math::sin_pi(x)/pi<Real>();
+    }
+
+    Real prime(Real t) const {
+        using boost::math::constants::pi;
+        using std::isfinite;
+        using std::floor;
+
+        Real x = (t - m_t0)/m_h;
+        if (ceil(x) == x) {
+            Real s = 0;
+            long j = static_cast<long>(x);
+            long n = m_y.size();
+            for (long i = 0; i < n; ++i)
+            {
+                if (j - i != 0)
+                {
+                    s += m_y[i]/(j-i);
+                }
+                // else derivative of sinc at zero is zero.
+            }
+            if (j & 1) {
+                s /= -m_h;
+            } else {
+                s /= m_h;
+            }
+            return s;
+        }
+        Real z = x;
+        auto it = m_y.begin();
+        Real cospix = boost::math::cos_pi(x);
+        Real sinpix_div_pi = boost::math::sin_pi(x)/pi<Real>();
+
+        Real s = 0;
+        auto end = m_y.end();
+        while(it != end)
+        {
+            s += (*it++)*(z*cospix - sinpix_div_pi)/(z*z);
+            z -= 1;
+        }
+
+        return s/m_h;
+    }
+
+
+
+    Real operator[](size_t i) const {
+        if (i & 1)
+        {
+            return -m_y[i];
+        }
+        return m_y[i];
+    }
+
+    RandomAccessContainer&& return_data() {
+        for (size_t i = 1; i < m_y.size(); i += 2)
+        {
+            m_y[i] = -m_y[i];
+        }
+        return std::move(m_y);
+    }
+
+
+private:
+    RandomAccessContainer m_y;
+    Real m_t0;
+    Real m_h;
+};
+}}}}
+#endif

boost/math/interpolators/makima.hpp

+// 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)
+
+// See: https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
+// And: https://doi.org/10.1145/321607.321609
+
+#ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
+#define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
+#include <memory>
+#include <cmath>
+#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template<class RandomAccessContainer>
+class makima {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    makima(RandomAccessContainer && x, RandomAccessContainer && y,
+           Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+           Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
+    {
+        using std::isnan;
+        using std::abs;
+        if (x.size() < 4)
+        {
+            throw std::domain_error("Must be at least four data points.");
+        }
+        RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
+        Real m2 = (y[3]-y[2])/(x[3]-x[2]);
+        Real m1 = (y[2]-y[1])/(x[2]-x[1]);
+        Real m0 = (y[1]-y[0])/(x[1]-x[0]);
+        // Quadratic extrapolation: m_{-1} = 2m_0 - m_1:
+        Real mm1 = 2*m0 - m1;
+        // Quadratic extrapolation: m_{-2} = 2*m_{-1}-m_0:
+        Real mm2 = 2*mm1 - m0;
+        Real w1 = abs(m1-m0) + abs(m1+m0)/2;
+        Real w2 = abs(mm1-mm2) + abs(mm1+mm2)/2;
+        if (isnan(left_endpoint_derivative))
+        {
+            s[0] = (w1*mm1 + w2*m0)/(w1+w2);
+            if (isnan(s[0]))
+            {
+                s[0] = 0;
+            }
+        }
+        else
+        {
+            s[0] = left_endpoint_derivative;
+        }
+
+        w1 = abs(m2-m1) + abs(m2+m1)/2;
+        w2 = abs(m0-mm1) + abs(m0+mm1)/2;
+        s[1] = (w1*m0 + w2*m1)/(w1+w2);
+        if (isnan(s[1])) {
+            s[1] = 0;
+        }
+
+        for (decltype(s.size()) i = 2; i < s.size()-2; ++i) {
+            Real mim2 = (y[i-1]-y[i-2])/(x[i-1]-x[i-2]);
+            Real mim1 = (y[i  ]-y[i-1])/(x[i  ]-x[i-1]);
+            Real mi   = (y[i+1]-y[i  ])/(x[i+1]-x[i  ]);
+            Real mip1 = (y[i+2]-y[i+1])/(x[i+2]-x[i+1]);
+            w1 = abs(mip1-mi) + abs(mip1+mi)/2;
+            w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
+            s[i] = (w1*mim1 + w2*mi)/(w1+w2);
+            if (isnan(s[i])) {
+                s[i] = 0;
+            }
+        }
+        // Quadratic extrapolation at the other end:
+        
+        decltype(s.size()) n = s.size();
+        Real mnm4 = (y[n-3]-y[n-4])/(x[n-3]-x[n-4]);
+        Real mnm3 = (y[n-2]-y[n-3])/(x[n-2]-x[n-3]);
+        Real mnm2 = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
+        Real mnm1 = 2*mnm2 - mnm3;
+        Real mn = 2*mnm1 - mnm2;
+        w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
+        w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
+
+        s[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
+        if (isnan(s[n-2])) {
+            s[n-2] = 0;
+        }
+
+        w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
+        w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
+
+
+        if (isnan(right_endpoint_derivative))
+        {
+            s[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
+            if (isnan(s[n-1])) {
+                s[n-1] = 0;
+            }
+        }
+        else
+        {
+            s[n-1] = right_endpoint_derivative;
+        }
+
+        impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
+    }
+
+    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 makima & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    void push_back(Real x, Real y) {
+        using std::abs;
+        using std::isnan;
+        if (x <= impl_->x_.back()) {
+             throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
+        }
+        impl_->x_.push_back(x);
+        impl_->y_.push_back(y);
+        impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
+        // dydx_[n-2] was computed by extrapolation. Now dydx_[n-2] -> dydx_[n-3], and it can be computed by the same formula.
+        decltype(impl_->size()) n = impl_->size();
+        auto i = n - 3;
+        Real mim2 = (impl_->y_[i-1]-impl_->y_[i-2])/(impl_->x_[i-1]-impl_->x_[i-2]);
+        Real mim1 = (impl_->y_[i  ]-impl_->y_[i-1])/(impl_->x_[i  ]-impl_->x_[i-1]);
+        Real mi   = (impl_->y_[i+1]-impl_->y_[i  ])/(impl_->x_[i+1]-impl_->x_[i  ]);
+        Real mip1 = (impl_->y_[i+2]-impl_->y_[i+1])/(impl_->x_[i+2]-impl_->x_[i+1]);
+        Real w1 = abs(mip1-mi) + abs(mip1+mi)/2;
+        Real w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
+        impl_->dydx_[i] = (w1*mim1 + w2*mi)/(w1+w2);
+        if (isnan(impl_->dydx_[i])) {
+            impl_->dydx_[i] = 0;
+        }
+
+        Real mnm4 = (impl_->y_[n-3]-impl_->y_[n-4])/(impl_->x_[n-3]-impl_->x_[n-4]);
+        Real mnm3 = (impl_->y_[n-2]-impl_->y_[n-3])/(impl_->x_[n-2]-impl_->x_[n-3]);
+        Real mnm2 = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1]-impl_->x_[n-2]);
+        Real mnm1 = 2*mnm2 - mnm3;
+        Real mn = 2*mnm1 - mnm2;
+        w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
+        w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
+
+        impl_->dydx_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
+        if (isnan(impl_->dydx_[n-2])) {
+            impl_->dydx_[n-2] = 0;
+        }
+
+        w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
+        w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
+
+        impl_->dydx_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
+        if (isnan(impl_->dydx_[n-1])) {
+            impl_->dydx_[n-1] = 0;
+        }
+    }
+
+private:
+    std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+}
+#endif
\ No newline at end of file

boost/math/interpolators/pchip.hpp

+// 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)
+
+#ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP
+#define BOOST_MATH_INTERPOLATORS_PCHIP_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template<class RandomAccessContainer>
+class pchip {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+
+    pchip(RandomAccessContainer && x, RandomAccessContainer && y,
+          Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
+          Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
+    {
+        if (x.size() < 4)
+        {
+            throw std::domain_error("Must be at least four data points.");
+        }
+        RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
+        if (isnan(left_endpoint_derivative))
+        {
+            // O(h) finite difference derivative:
+            // This, I believe, is the only derivative guaranteed to be monotonic:
+            s[0] = (y[1]-y[0])/(x[1]-x[0]);
+        }
+        else
+        {
+            s[0] = left_endpoint_derivative;
+        }
+
+        for (decltype(s.size()) k = 1; k < s.size()-1; ++k) {
+            Real hkm1 = x[k] - x[k-1];
+            Real dkm1 = (y[k] - y[k-1])/hkm1;
+
+            Real hk = x[k+1] - x[k];
+            Real dk = (y[k+1] - y[k])/hk;
+            Real w1 = 2*hk + hkm1;
+            Real w2 = hk + 2*hkm1;
+            if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
+            {
+                s[k] = 0;
+            }
+            else
+            {
+                s[k] = (w1+w2)/(w1/dkm1 + w2/dk);
+            }
+
+        }
+        // Quadratic extrapolation at the other end:
+        auto n = s.size();
+        if (isnan(right_endpoint_derivative))
+        {
+            s[n-1] = (y[n-1]-y[n-2])/(x[n-1] - x[n-2]);
+        }
+        else
+        {
+            s[n-1] = right_endpoint_derivative;
+        }
+        impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
+    }
+
+    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 pchip & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    void push_back(Real x, Real y) {
+        using std::abs;
+        using std::isnan;
+        if (x <= impl_->x_.back()) {
+             throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
+        }
+        impl_->x_.push_back(x);
+        impl_->y_.push_back(y);
+        impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
+        auto n = impl_->size();
+        impl_->dydx_[n-1] = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1] - impl_->x_[n-2]);
+        // Now fix s_[n-2]:
+        auto k = n-2;
+        Real hkm1 = impl_->x_[k] - impl_->x_[k-1];
+        Real dkm1 = (impl_->y_[k] - impl_->y_[k-1])/hkm1;
+
+        Real hk = impl_->x_[k+1] - impl_->x_[k];
+        Real dk = (impl_->y_[k+1] - impl_->y_[k])/hk;
+        Real w1 = 2*hk + hkm1;
+        Real w2 = hk + 2*hkm1;
+        if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
+        {
+            impl_->dydx_[k] = 0;
+        }
+        else
+        {
+            impl_->dydx_[k] = (w1+w2)/(w1/dkm1 + w2/dk);
+        }
+    }
+
+private:
+    std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+}
+#endif
\ No newline at end of file

boost/math/interpolators/quintic_hermite.hpp

+/*
+ * 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)
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP
+#define BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP
+#include <algorithm>
+#include <stdexcept>
+#include <memory>
+#include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template<class RandomAccessContainer>
+class quintic_hermite {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2)
+     : impl_(std::make_shared<detail::quintic_hermite_detail<RandomAccessContainer>>(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);
+    }
+
+    Real double_prime(Real x) const
+    {
+        return impl_->double_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);
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+private:
+    std::shared_ptr<detail::quintic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+template<class RandomAccessContainer>
+class cardinal_quintic_hermite {
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    cardinal_quintic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, Real x0, Real dx)
+     : impl_(std::make_shared<detail::cardinal_quintic_hermite_detail<RandomAccessContainer>>(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx))
+    {}
+
+    inline Real operator()(Real x) const {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const {
+        return impl_->prime(x);
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        return impl_->double_prime(x);
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+private:
+    std::shared_ptr<detail::cardinal_quintic_hermite_detail<RandomAccessContainer>> impl_;
+};
+
+template<class RandomAccessContainer>
+class cardinal_quintic_hermite_aos {
+public:
+    using Point = typename RandomAccessContainer::value_type;
+    using Real = typename Point::value_type;
+    cardinal_quintic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx)
+     : impl_(std::make_shared<detail::cardinal_quintic_hermite_detail_aos<RandomAccessContainer>>(std::move(data), x0, dx))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        return impl_->double_prime(x);
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+private:
+    std::shared_ptr<detail::cardinal_quintic_hermite_detail_aos<RandomAccessContainer>> impl_;
+};
+
+
+}
+#endif

boost/math/interpolators/septic_hermite.hpp

+/*
+ * 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)
+ */
+#ifndef BOOST_MATH_INTERPOLATORS_SEPTIC_HERMITE_HPP
+#define BOOST_MATH_INTERPOLATORS_SEPTIC_HERMITE_HPP
+#include <algorithm>
+#include <stdexcept>
+#include <memory>
+#include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
+
+namespace boost::math::interpolators {
+
+template<class RandomAccessContainer>
+class septic_hermite
+{
+public:
+    using Real = typename RandomAccessContainer::value_type;
+    septic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, 
+                   RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3)
+     : impl_(std::make_shared<detail::septic_hermite_detail<RandomAccessContainer>>(std::move(x), 
+     std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3)))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        return impl_->double_prime(x);
+    }
+
+    friend std::ostream& operator<<(std::ostream & os, const septic_hermite & m)
+    {
+        os << *m.impl_;
+        return os;
+    }
+
+    int64_t bytes() const
+    {
+        return impl_->bytes() + sizeof(impl_);
+    }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
+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))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    inline Real double_prime(Real x) const
+    {
+        return impl_->double_prime(x);
+    }
+
+    int64_t bytes() const
+    {
+        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_;
+};
+
+
+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))
+    {}
+
+    inline Real operator()(Real x) const
+    {
+        return impl_->operator()(x);
+    }
+
+    inline Real prime(Real x) const
+    {
+        return impl_->prime(x);
+    }
+
+    inline Real double_prime(Real x) const 
+    {
+        return impl_->double_prime(x);
+    }
+
+    int64_t bytes() const
+    {
+        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_;
+};
+
+}
+#endif 
\ No newline at end of file

boost/math/interpolators/vector_barycentric_rational.hpp

+/*
+ *  Copyright Nick Thompson, 2019
+ *  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)
+ *
+ *  Exactly the same as barycentric_rational.hpp, but delivers values in $\mathbb{R}^n$.
+ *  In some sense this is trivial, since each component of the vector is computed in exactly the same
+ *  as would be computed by barycentric_rational.hpp. But this is a bit more efficient and convenient.
+ */
+
+#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP
+#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP
+
+#include <memory>
+#include <boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp>
+
+namespace boost{ namespace math{
+
+template<class TimeContainer, class SpaceContainer>
+class vector_barycentric_rational
+{
+public:
+    using Real = typename TimeContainer::value_type;
+    using Point = typename SpaceContainer::value_type;
+    vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order = 3);
+
+    void operator()(Point& x, Real t) const;
+
+    // I have validated using google benchmark that returning a value is no more expensive populating it,
+    // at least for Eigen vectors with known size at compile-time.
+    // This is kinda a weird thing to discover since it goes against the advice of basically every high-performance computing book.
+    Point operator()(Real t) const {
+        Point p;
+        this->operator()(p, t);
+        return p;
+    }
+
+    void prime(Point& dxdt, Real t) const {
+        Point x;
+        m_imp->eval_with_prime(x, dxdt, t);
+    }
+
+    Point prime(Real t) const {
+        Point p;
+        this->prime(p, t);
+        return p;
+    }
+
+    void eval_with_prime(Point& x, Point& dxdt, Real t) const {
+        m_imp->eval_with_prime(x, dxdt, t);
+        return;
+    }
+
+    std::pair<Point, Point> eval_with_prime(Real t) const {
+        Point x;
+        Point dxdt;
+        m_imp->eval_with_prime(x, dxdt, t);
+        return {x, dxdt};
+    }
+
+private:
+    std::shared_ptr<detail::vector_barycentric_rational_imp<TimeContainer, SpaceContainer>> m_imp;
+};
+
+
+template <class TimeContainer, class SpaceContainer>
+vector_barycentric_rational<TimeContainer, SpaceContainer>::vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order):
+ m_imp(std::make_shared<detail::vector_barycentric_rational_imp<TimeContainer, SpaceContainer>>(std::move(times), std::move(points), approximation_order))
+{
+    return;
+}
+
+template <class TimeContainer, class SpaceContainer>
+void vector_barycentric_rational<TimeContainer, SpaceContainer>::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const
+{
+    m_imp->operator()(p, t);
+    return;
+}
+
+}}
+#endif

boost/math/interpolators/whittaker_shannon.hpp

+// Copyright Nick Thompson, 2019
+// 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)
+#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
+#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
+#include <memory>
+#include <boost/math/interpolators/detail/whittaker_shannon_detail.hpp>
+
+namespace boost { namespace math { namespace interpolators {
+
+template<class RandomAccessContainer>
+class whittaker_shannon {
+public:
+
+    using Real = typename RandomAccessContainer::value_type;
+    whittaker_shannon(RandomAccessContainer&& y, Real const & t0, Real const & h)
+     : m_impl(std::make_shared<detail::whittaker_shannon_detail<RandomAccessContainer>>(std::move(y), t0, h))
+    {}
+
+    inline Real operator()(Real t) const
+    {
+        return m_impl->operator()(t);
+    }
+
+    inline Real prime(Real t) const
+    {
+        return m_impl->prime(t);
+    }
+
+    inline Real operator[](size_t i) const
+    {
+        return m_impl->operator[](i);
+    }
+
+    RandomAccessContainer&& return_data()
+    {
+        return m_impl->return_data();
+    }
+
+
+private:
+    std::shared_ptr<detail::whittaker_shannon_detail<RandomAccessContainer>> m_impl;
+};
+}}}
+#endif