diff --git a/.travis.yml b/.travis.yml
index e61d94154..5d46c76a2 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -571,7 +571,7 @@ matrix:
             - clang-6.0
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-trusty-6.0
+            - llvm-toolchain-xenial-6.0
 
     - os: linux
       compiler: clang++-6.0
@@ -582,7 +582,7 @@ matrix:
             - clang-6.0
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-trusty-6.0
+            - llvm-toolchain-xenial-6.0
 
     - os: linux
       compiler: clang++-6.0
@@ -593,7 +593,7 @@ matrix:
             - clang-6.0
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-trusty-6.0
+            - llvm-toolchain-xenial-6.0
 
     - os: linux
       compiler: clang++-6.0
@@ -604,7 +604,7 @@ matrix:
             - clang-6.0
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-trusty-6.0
+            - llvm-toolchain-xenial-6.0
 
     - os: linux
       compiler: clang++-6.0
@@ -615,7 +615,7 @@ matrix:
             - clang-6.0
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-trusty-6.0
+            - llvm-toolchain-xenial-6.0
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun
diff --git a/doc/interpolators/vector_barycentric_rational.qbk b/doc/interpolators/vector_barycentric_rational.qbk
new file mode 100644
index 000000000..335ecf372
--- /dev/null
+++ b/doc/interpolators/vector_barycentric_rational.qbk
@@ -0,0 +1,76 @@
+[/
+  Copyright 2019 Nick Thompson
+
+  Distributed under 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).
+]
+
+[section:vector_barycentric Vector-valued Barycentric Rational Interpolation]
+
+[heading Synopsis]
+
+``
+#include <boost/math/interpolators/vector_barycentric_rational.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;
+
+    Point operator()(Real t) const;
+
+    void prime(Point& dxdt, Real t) const;
+
+    Point prime(Real t);
+
+    void eval_with_prime(Point& x, Point& dxdt, Real t) const;
+
+    std::pair<Point, Point> eval_with_prime(Real t) const;
+};
+
+}}
+``
+
+[heading Description]
+
+The /n/ dimensional vector-valued barycentric rational interpolator is exactly the same as /n/ scalar-valued barycentric rational interpolators.
+This is provided primarily for convenience and a slight improvement in efficiency over using /n/ different rational interpolators and combining their results.
+
+Use of the class requires a `Point`-type which has size known at compile time.
+These requirements are satisfied by (for example) `Eigen::Vector2d`s and `std::array<Real, N>` classes.
+The call to the constructor computes the weights:
+
+    using boost::math::vector_barycentric_rational;
+    std::vector<double> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    // initialize t and y . . .
+    vector_barycentric_rational<decltype(t), decltype(y)> interpolant(std::move(t), std::move(y));
+
+To evaluate the interpolant, use
+
+    double t = 2.3;
+    Eigen::Vector2d y = interpolant(t);
+
+If you want to populate a vector passed into the interpolant, rather than get it returned, that syntax is supported:
+
+    Eigen::Vector2d y;
+    interpolant(y, t);
+
+We tested this with `Eigen::Vector`s and found no performance benefit, but other `Point`-types might not be the same.
+
+To evaluate the derivative of the interpolant use
+
+    auto [y, y_prime] = interpolant.eval_with_prime(x);
+
+Computation of the derivative requires evaluation, so if you can try to use both values at once.
+
+
+[endsect] [/section:vector_barycentric Vector Barycentric Rational Interpolation]
diff --git a/doc/math.qbk b/doc/math.qbk
index e64256403..a3ce216ad 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -662,6 +662,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [mathpart interpolation Interpolation]
 [include interpolators/cubic_b_spline.qbk]
 [include interpolators/barycentric_rational_interpolation.qbk]
+[include interpolators/vector_barycentric_rational.qbk]
 [include interpolators/catmull_rom.qbk]
 [endmathpart]
 
diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp
new file mode 100644
index 000000000..2f8d177a3
--- /dev/null
+++ b/include/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/include/boost/math/interpolators/vector_barycentric_rational.hpp b/include/boost/math/interpolators/vector_barycentric_rational.hpp
new file mode 100644
index 000000000..8289799bc
--- /dev/null
+++ b/include/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/test/Jamfile.v2 b/test/Jamfile.v2
index ab67fe9bb..c16598f7e 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -899,6 +899,7 @@ test-suite distribution_tests :
 test-suite misc :
    [ run test_print_info_on_type.cpp  ]
    [ run test_barycentric_rational.cpp ../../test/build//boost_unit_test_framework : : :  [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ]  ]
+   [ run test_vector_barycentric_rational.cpp ../../test/build//boost_unit_test_framework : : :  [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ]  [ check-target-builds ../../multiprecision/config//has_eigen : : <build>no ] ]
    [ run test_constant_generate.cpp : : : release <define>USE_CPP_FLOAT=1 <exception-handling>off:<build>no  ]
    [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : :  [ requires cxx11_smart_ptr cxx11_defaulted_functions ] <debug-symbols>off <toolset>msvc:<cxxflags>/bigobj release ]
    [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ]
diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp
new file mode 100644
index 000000000..6250fa897
--- /dev/null
+++ b/test/test_vector_barycentric_rational.cpp
@@ -0,0 +1,385 @@
+// 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)
+
+#define BOOST_TEST_MODULE vector_barycentric_rational
+
+#include <cmath>
+#include <random>
+#include <array>
+#include <Eigen/Dense>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
+#include <boost/type_index.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/math/interpolators/barycentric_rational.hpp>
+#include <boost/math/interpolators/vector_barycentric_rational.hpp>
+
+using std::sqrt;
+using std::abs;
+using std::numeric_limits;
+
+template<class Real>
+void test_agreement_with_1d()
+{
+    std::cout << "Testing with 1D interpolation on type "
+              << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+    std::mt19937 gen(4723);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    std::vector<Real> t_copy0 = t;
+    std::vector<Real> t_copy1 = t;
+
+    std::vector<Real> y_copy0(y.size());
+    std::vector<Real> y_copy1(y.size());
+    for (size_t i = 0; i < y.size(); ++i) {
+        y_copy0[i] = y[i][0];
+        y_copy1[i] = y[i][1];
+    }
+
+    boost::random::uniform_real_distribution<Real> dis2(t[0], t[t.size()-1]);
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+    boost::math::barycentric_rational<Real> scalar_interpolator0(std::move(t_copy0), std::move(y_copy0));
+    boost::math::barycentric_rational<Real> scalar_interpolator1(std::move(t_copy1), std::move(y_copy1));
+
+
+    Eigen::Vector2d z;
+
+    size_t samples = 0;
+    while (samples++ < 1000)
+    {
+        Real t = dis2(gen);
+        interpolator(z, t);
+        BOOST_CHECK_CLOSE(z[0], scalar_interpolator0(t), 2*numeric_limits<Real>::epsilon());
+        BOOST_CHECK_CLOSE(z[1], scalar_interpolator1(t), 2*numeric_limits<Real>::epsilon());
+    }
+}
+
+
+template<class Real>
+void test_interpolation_condition_eigen()
+{
+    std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type "
+              << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+    std::mt19937 gen(4723);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+
+    Eigen::Vector2d z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        interpolator(z, t_copy[i]);
+        BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
+        BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
+    }
+}
+
+template<class Real>
+void test_interpolation_condition_std_array()
+{
+    std::cout << "Testing interpolation condition for barycentric interpolation on std::array vectors of type "
+              << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+    std::mt19937 gen(4723);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<std::array<Real, 2>> y(100);
+    t[0] = dis(gen);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<std::array<Real, 2>> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+
+    std::array<Real, 2> z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        interpolator(z, t_copy[i]);
+        BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
+        BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
+    }
+}
+
+
+template<class Real>
+void test_interpolation_condition_ublas()
+{
+    std::cout << "Testing interpolation condition for barycentric interpolation ublas vectors of type "
+              << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+    std::mt19937 gen(4723);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<boost::numeric::ublas::vector<Real>> y(100);
+    t[0] = dis(gen);
+    y[0].resize(2);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i].resize(2);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<Real> t_copy = t;
+    std::vector<boost::numeric::ublas::vector<Real>> y_copy = y;
+
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+
+    boost::numeric::ublas::vector<Real> z(2);
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        interpolator(z, t_copy[i]);
+        BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
+        BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
+    }
+}
+
+template<class Real>
+void test_interpolation_condition_high_order()
+{
+    std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+    std::mt19937 gen(5);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
+
+    Eigen::Vector2d z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        interpolator(z, t_copy[i]);
+        BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
+        BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
+    }
+}
+
+
+template<class Real>
+void test_constant_eigen()
+{
+    std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on Eigen vectors of type "
+              << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+
+    std::mt19937 gen(6);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    Real constant0 = dis(gen);
+    Real constant1 = dis(gen);
+    y[0][0] = constant0;
+    y[0][1] = constant1;
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = constant0;
+        y[i][1] = constant1;
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+
+    Eigen::Vector2d z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
+        Real t = t_copy[i] + dis(gen);
+        z = interpolator(t);
+        BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
+        Eigen::Vector2d zprime = interpolator.prime(t);
+        Real zero_0 = zprime[0];
+        Real zero_1 = zprime[1];
+        BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
+    }
+}
+
+
+template<class Real>
+void test_constant_std_array()
+{
+    std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on std::array vectors of type "
+              << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+
+    std::mt19937 gen(6);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<std::array<Real, 2>> y(100);
+    t[0] = dis(gen);
+    Real constant0 = dis(gen);
+    Real constant1 = dis(gen);
+    y[0][0] = constant0;
+    y[0][1] = constant1;
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = constant0;
+        y[i][1] = constant1;
+    }
+
+    std::vector<std::array<Real,2>> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
+
+    std::array<Real, 2> z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
+        Real t = t_copy[i] + dis(gen);
+        z = interpolator(t);
+        BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
+        std::array<Real, 2> zprime = interpolator.prime(t);
+        Real zero_0 = zprime[0];
+        Real zero_1 = zprime[1];
+        BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
+    }
+}
+
+
+template<class Real>
+void test_constant_high_order()
+{
+    std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+
+    std::mt19937 gen(6);
+    boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    Real constant0 = dis(gen);
+    Real constant1 = dis(gen);
+    y[0][0] = constant0;
+    y[0][1] = constant1;
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = constant0;
+        y[i][1] = constant1;
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
+
+    Eigen::Vector2d z;
+    for (size_t i = 0; i < t_copy.size(); ++i)
+    {
+        // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
+        Real t = t_copy[i] + dis(gen);
+        z = interpolator(t);
+        BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
+        Eigen::Vector2d zprime = interpolator.prime(t);
+        Real zero_0 = zprime[0];
+        Real zero_1 = zprime[1];
+        BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
+        BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
+    }
+}
+
+
+template<class Real>
+void test_weights()
+{
+    std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+
+    std::mt19937 gen(9);
+    boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
+    std::vector<Real> t(100);
+    std::vector<Eigen::Vector2d> y(100);
+    t[0] = dis(gen);
+    y[0][0] = dis(gen);
+    y[0][1] = dis(gen);
+    for (size_t i = 1; i < t.size(); ++i)
+    {
+        t[i] = t[i-1] + dis(gen);
+        y[i][0] = dis(gen);
+        y[i][1] = dis(gen);
+    }
+
+    std::vector<Eigen::Vector2d> y_copy = y;
+    std::vector<Real> t_copy = t;
+    boost::math::detail::vector_barycentric_rational_imp<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 1);
+
+    for (size_t i = 1; i < t_copy.size() - 1; ++i)
+    {
+        Real w = interpolator.weight(i);
+        Real w_expect = 1/(t_copy[i] - t_copy[i - 1]) + 1/(t_copy[i+1] - t_copy[i]);
+        if (i % 2 == 0)
+        {
+            BOOST_CHECK_CLOSE(w, -w_expect, 0.00001);
+        }
+        else
+        {
+            BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(vector_barycentric_rational)
+{
+    test_weights<double>();
+    test_constant_eigen<double>();
+    test_constant_std_array<double>();
+    test_constant_high_order<double>();
+    test_interpolation_condition_eigen<double>();
+    test_interpolation_condition_ublas<double>();
+    test_interpolation_condition_std_array<double>();
+    test_interpolation_condition_high_order<double>();
+    test_agreement_with_1d<double>();
+}