ceph/src/boost/boost/math/interpolators/cubic_b_spline.hpp

   1 // Copyright Nick Thompson, 2017
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 // This implements the compactly supported cubic b spline algorithm described in
   8 // Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
   9 // Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
  10
  11 // Let f be the function we are trying to interpolate, and s be the interpolating spline.
  12 // The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
  13 // The order of accuracy depends on the regularity of the f, however, assuming f is
  14 // four-times continuously differentiable, the error is of O(h^4).
  15 // In addition, we can differentiate the spline and obtain a good interpolant for f'.
  16 // The main restriction of this method is that the samples of f must be evenly spaced.
  17 // Look for barycentric rational interpolation for non-evenly sampled data.
  18 // Properties:
  19 // - s(x_j) = f(x_j)
  20 // - All cubic polynomials interpolated exactly
  21
  22 #ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
  23 #define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
  24
  25 #include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
  26 #include <boost/config/header_deprecated.hpp>
  27
  28 BOOST_HEADER_DEPRECATED("<boost/math/interpolators/cardinal_cubic_b_spline.hpp>");
  29
  30 namespace boost{ namespace math{
  31
  32 template <class Real>
  33 class cubic_b_spline
  34 {
  35 public:
  36     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
  37     // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
  38     template <class BidiIterator>
  39     cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
  40                    Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
  41                    Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
  42     cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
  43        Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
  44        Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
  45
  46     cubic_b_spline() = default;
  47     Real operator()(Real x) const;
  48
  49     Real prime(Real x) const;
  50
  51     Real double_prime(Real x) const;
  52
  53 private:
  54     std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
  55 };
  56
  57 template<class Real>
  58 cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
  59                                      Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
  60 {
  61 }
  62
  63 template <class Real>
  64 template <class BidiIterator>
  65 cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
  66    Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
  67 {
  68 }
  69
  70 template<class Real>
  71 Real cubic_b_spline<Real>::operator()(Real x) const
  72 {
  73     return m_imp->operator()(x);
  74 }
  75
  76 template<class Real>
  77 Real cubic_b_spline<Real>::prime(Real x) const
  78 {
  79     return m_imp->prime(x);
  80 }
  81
  82 template<class Real>
  83 Real cubic_b_spline<Real>::double_prime(Real x) const
  84 {
  85     return m_imp->double_prime(x);
  86 }
  87
  88
  89 }}
  90 #endif