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