boost/math/interpolators/barycentric_rational.hpp

   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
  49     std::vector<Real>&& return_x()
  50     {
  51         return m_imp->return_x();
  52     }
  53
  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>
  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>
  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>
  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>
  86 Real barycentric_rational<Real>::operator()(Real x) const
  87 {
  88     return m_imp->operator()(x);
  89 }
  90
  91 template<class Real>
  92 Real barycentric_rational<Real>::prime(Real x) const
  93 {
  94     return m_imp->prime(x);
  95 }
  96
  97
  98 }}
  99 #endif