1 /*
2 * Copyright Nick Thompson, 2017
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 *
7 * Given N samples (t_i, y_i) which are irregularly spaced, this routine constructs an
8 * interpolant s which is constructed in O(N) time, occupies O(N) space, and can be evaluated in O(N) time.
9 * The interpolation is stable, unless one point is incredibly close to another, and the next point is incredibly far.
10 * The measure of this stability is the "local mesh ratio", which can be queried from the routine.
11 * Pictorially, the following t_i spacing is bad (has a high local mesh ratio)
12 * || | | | |
13 * and this t_i spacing is good (has a low local mesh ratio)
14 * | | | | | | | | | |
15 *
16 *
17 * If f is C^{d+2}, then the interpolant is O(h^(d+1)) accurate, where d is the interpolation order.
18 * A disadvantage of this interpolant is that it does not reproduce rational functions; for example, 1/(1+x^2) is not interpolated exactly.
19 *
20 * References:
21 * Floater, Michael S., and Kai Hormann. "Barycentric rational interpolation with no poles and high rates of approximation."
22 * Numerische Mathematik 107.2 (2007): 315-331.
23 * Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
24 */
25
26 #ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
27 #define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
28
29 #include <memory>
30 #include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
31
32 namespace boost{ namespace math{
33
34 template<class Real>
35 class barycentric_rational
36 {
37 public:
38 barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
39
40 barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
41
42 template <class InputIterator1, class InputIterator2>
43 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);
44
45 Real operator()(Real x) const;
46
47 Real prime(Real x) const;
48
return_x()49 std::vector<Real>&& return_x()
50 {
51 return m_imp->return_x();
52 }
53
return_y()54 std::vector<Real>&& return_y()
55 {
56 return m_imp->return_y();
57 }
58
59 private:
60 std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
61 };
62
63 template <class Real>
barycentric_rational(const Real * const x,const Real * const y,size_t n,size_t approximation_order)64 barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
65 m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
66 {
67 return;
68 }
69
70 template <class Real>
barycentric_rational(std::vector<Real> && x,std::vector<Real> && y,size_t approximation_order)71 barycentric_rational<Real>::barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order):
72 m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(std::move(x), std::move(y), approximation_order))
73 {
74 return;
75 }
76
77
78 template <class Real>
79 template <class InputIterator1, class InputIterator2>
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 *)80 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*)
81 : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
82 {
83 }
84
85 template<class Real>
operator ()(Real x) const86 Real barycentric_rational<Real>::operator()(Real x) const
87 {
88 return m_imp->operator()(x);
89 }
90
91 template<class Real>
prime(Real x) const92 Real barycentric_rational<Real>::prime(Real x) const
93 {
94 return m_imp->prime(x);
95 }
96
97
98 }}
99 #endif
100