Imported Upstream version 1.71.0upstream/1.71.0

7 files changed, 753 insertions, 22 deletions

diff --git a/boost/math/interpolators/cardinal_quadratic_b_spline.hpp b/boost/math/interpolators/cardinal_quadratic_b_spline.hpp
new file mode 100644
index 0000000000..a5b150f2f1
--- /dev/null
+++ b/boost/math/interpolators/cardinal_quadratic_b_spline.hpp
@@ -0,0 +1,57 @@
+// 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
diff --git a/boost/math/interpolators/catmull_rom.hpp b/boost/math/interpolators/catmull_rom.hpp
index 05ad666de2..20207871db 100644
--- a/boost/math/interpolators/catmull_rom.hpp
+++ b/boost/math/interpolators/catmull_rom.hpp
@@ -13,6 +13,25 @@
 #include <vector>
 #include <algorithm>
 #include <iterator>
+#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{
 
@@ -22,7 +41,7 @@ namespace boost{ namespace math{
         typename Point::value_type alpha_distance(Point const & p1, Point const & p2, typename Point::value_type alpha)
         {
             using std::pow;
-            using std::size;
+            using std_workaround::size;
             typename Point::value_type dsq = 0;
             for (size_t i = 0; i < size(p1); ++i)
             {
@@ -33,44 +52,45 @@ namespace boost{ namespace math{
         }
     }
 
-template <class Point>
+template <class Point, class RandomAccessContainer = std::vector<Point> >
 class catmull_rom
 {
+   typedef typename Point::value_type value_type;
 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(RandomAccessContainer&& points, bool closed = false, value_type alpha = (value_type) 1/ (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) {}
+    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) {}
 
-    typename Point::value_type max_parameter() const
+    value_type max_parameter() const
     {
         return m_max_s;
     }
 
-    typename Point::value_type parameter_at_point(size_t i) const
+    value_type parameter_at_point(size_t i) const
     {
         return m_s[i+1];
     }
 
-    Point operator()(const typename Point::value_type s) const;
+    Point operator()(const value_type s) const;
 
-    Point prime(const typename Point::value_type s) const;
+    Point prime(const value_type s) const;
 
-    std::vector<Point>&& get_points()
+    RandomAccessContainer&& 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;
+    RandomAccessContainer m_pnts;
+    std::vector<value_type> m_s;
+    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))
+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))
 {
-    size_t num_pnts = m_pnts.size();
+    std::size_t num_pnts = m_pnts.size();
     //std::cout << "Number of points = " << num_pnts << "\n";
     if (num_pnts < 4)
     {
@@ -90,7 +110,7 @@ catmull_rom<Point>::catmull_rom(std::vector<Point>&& points, bool closed, typena
     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)
+    for (std::ptrdiff_t i = num_pnts-1; i >= 0; --i)
     {
         m_pnts[i+1] = m_pnts[i];
     }
@@ -122,10 +142,10 @@ catmull_rom<Point>::catmull_rom(std::vector<Point>&& points, bool closed, typena
 }
 
 
-template<class Point>
-Point catmull_rom<Point>::operator()(const typename Point::value_type s) const
+template<class Point, class RandomAccessContainer >
+Point catmull_rom<Point, RandomAccessContainer>::operator()(const typename Point::value_type s) const
 {
-    using std::size;
+    using std_workaround::size;
     if (s < 0 || s > m_max_s)
     {
         throw std::domain_error("Parameter outside bounds.");
@@ -182,10 +202,10 @@ Point catmull_rom<Point>::operator()(const typename Point::value_type s) const
     return B1_or_C;
 }
 
-template<class Point>
-Point catmull_rom<Point>::prime(const typename Point::value_type s) const
+template<class Point, class RandomAccessContainer >
+Point catmull_rom<Point, RandomAccessContainer>::prime(const typename Point::value_type s) const
 {
-    using std::size;
+    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)
diff --git a/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp b/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp
new file mode 100644
index 0000000000..044ac72179
--- /dev/null
+++ b/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp
@@ -0,0 +1,206 @@
+// 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>
+
+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, γ be defined via t = t0 + kh + γh.
+        // Now find all j: |k-j+1+γ|< 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, γ be defined via t = t0 + kh + γh.
+        // Now find all j: |k-j+1+γ|< 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
diff --git a/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp
new file mode 100644
index 0000000000..2f8d177a35
--- /dev/null
+++ b/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp
@@ -0,0 +1,193 @@
+/*
+ *  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 <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 Real = typename TimeContainer::value_type;
+    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<Real>::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
diff --git a/boost/math/interpolators/detail/whittaker_shannon_detail.hpp b/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
new file mode 100644
index 0000000000..3d4ace397e
--- /dev/null
+++ b/boost/math/interpolators/detail/whittaker_shannon_detail.hpp
@@ -0,0 +1,126 @@
+// 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
diff --git a/boost/math/interpolators/vector_barycentric_rational.hpp b/boost/math/interpolators/vector_barycentric_rational.hpp
new file mode 100644
index 0000000000..8289799bc9
--- /dev/null
+++ b/boost/math/interpolators/vector_barycentric_rational.hpp
@@ -0,0 +1,82 @@
+/*
+ *  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
diff --git a/boost/math/interpolators/whittaker_shannon.hpp b/boost/math/interpolators/whittaker_shannon.hpp
new file mode 100644
index 0000000000..a5c5367a9d
--- /dev/null
+++ b/boost/math/interpolators/whittaker_shannon.hpp
@@ -0,0 +1,47 @@
+// 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