barycentric_rational.hpp 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899
  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." Numerische Mathematik 107.2 (2007): 315-331.
  22. * Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
  23. */
  24. #ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
  25. #define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_HPP
  26. #include <memory>
  27. #include <boost/math/interpolators/detail/barycentric_rational_detail.hpp>
  28. namespace boost{ namespace math{
  29. template<class Real>
  30. class barycentric_rational
  31. {
  32. public:
  33. barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order = 3);
  34. barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
  35. template <class InputIterator1, class InputIterator2>
  36. 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);
  37. Real operator()(Real x) const;
  38. Real prime(Real x) const;
  39. std::vector<Real>&& return_x()
  40. {
  41. return m_imp->return_x();
  42. }
  43. std::vector<Real>&& return_y()
  44. {
  45. return m_imp->return_y();
  46. }
  47. private:
  48. std::shared_ptr<detail::barycentric_rational_imp<Real>> m_imp;
  49. };
  50. template <class Real>
  51. barycentric_rational<Real>::barycentric_rational(const Real* const x, const Real* const y, size_t n, size_t approximation_order):
  52. m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(x, x + n, y, approximation_order))
  53. {
  54. return;
  55. }
  56. template <class Real>
  57. barycentric_rational<Real>::barycentric_rational(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order):
  58. m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(std::move(x), std::move(y), approximation_order))
  59. {
  60. return;
  61. }
  62. template <class Real>
  63. template <class InputIterator1, class InputIterator2>
  64. 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*)
  65. : m_imp(std::make_shared<detail::barycentric_rational_imp<Real>>(start_x, end_x, start_y, approximation_order))
  66. {
  67. }
  68. template<class Real>
  69. Real barycentric_rational<Real>::operator()(Real x) const
  70. {
  71. return m_imp->operator()(x);
  72. }
  73. template<class Real>
  74. Real barycentric_rational<Real>::prime(Real x) const
  75. {
  76. return m_imp->prime(x);
  77. }
  78. }}
  79. #endif