ceph/src/boost/boost/math/interpolators/detail/barycentric_rational_detail.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
   8 #ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
   9 #define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
  10
  11 #include <vector>
  12 #include <utility> // for std::move
  13 #include <algorithm> // for std::is_sorted
  14 #include <string>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/math/tools/assert.hpp>
  17
  18 namespace boost{ namespace math{ namespace interpolators { namespace detail{
  19
  20 template<class Real>
  21 class barycentric_rational_imp
  22 {
  23 public:
  24     template <class InputIterator1, class InputIterator2>
  25     barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
  26
  27     barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
  28
  29     Real operator()(Real x) const;
  30
  31     Real prime(Real x) const;
  32
  33     // The barycentric weights are not really that interesting; except to the unit tests!
  34     Real weight(size_t i) const { return m_w[i]; }
  35
  36     std::vector<Real>&& return_x()
  37     {
  38         return std::move(m_x);
  39     }
  40
  41     std::vector<Real>&& return_y()
  42     {
  43         return std::move(m_y);
  44     }
  45
  46 private:
  47
  48     void calculate_weights(size_t approximation_order);
  49
  50     std::vector<Real> m_x;
  51     std::vector<Real> m_y;
  52     std::vector<Real> m_w;
  53 };
  54
  55 template <class Real>
  56 template <class InputIterator1, class InputIterator2>
  57 barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
  58 {
  59     std::ptrdiff_t n = std::distance(start_x, end_x);
  60
  61     if (approximation_order >= (std::size_t)n)
  62     {
  63         throw std::domain_error("Approximation order must be < data length.");
  64     }
  65
  66     // Big sad memcpy.
  67     m_x.resize(n);
  68     m_y.resize(n);
  69     for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
  70     {
  71         // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
  72         if(boost::math::isnan(*start_x))
  73         {
  74             std::string msg = std::string("x[") + std::to_string(i) + "] is a NAN";
  75             throw std::domain_error(msg);
  76         }
  77
  78         if(boost::math::isnan(*start_y))
  79         {
  80            std::string msg = std::string("y[") + std::to_string(i) + "] is a NAN";
  81            throw std::domain_error(msg);
  82         }
  83
  84         m_x[i] = *start_x;
  85         m_y[i] = *start_y;
  86     }
  87     calculate_weights(approximation_order);
  88 }
  89
  90 template <class Real>
  91 barycentric_rational_imp<Real>::barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y))
  92 {
  93     BOOST_MATH_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates.");
  94     BOOST_MATH_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length.");
  95     BOOST_MATH_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1].");
  96     calculate_weights(approximation_order);
  97 }
  98
  99 template<class Real>
 100 void barycentric_rational_imp<Real>::calculate_weights(size_t approximation_order)
 101 {
 102     using std::abs;
 103     int64_t n = m_x.size();
 104     m_w.resize(n, 0);
 105     for(int64_t k = 0; k < n; ++k)
 106     {
 107         int64_t i_min = (std::max)(k - static_cast<int64_t>(approximation_order), static_cast<int64_t>(0));
 108         int64_t i_max = k;
 109         if (k >= n - (std::ptrdiff_t)approximation_order)
 110         {
 111             i_max = n - approximation_order - 1;
 112         }
 113
 114         for(int64_t i = i_min; i <= i_max; ++i)
 115         {
 116             Real inv_product = 1;
 117             int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
 118             for(int64_t j = i; j <= j_max; ++j)
 119             {
 120                 if (j == k)
 121                 {
 122                     continue;
 123                 }
 124
 125                 Real diff = m_x[k] - m_x[j];
 126                 using std::numeric_limits;
 127                 if (abs(diff) < (numeric_limits<Real>::min)())
 128                 {
 129                    std::string msg = std::string("Spacing between  x[")
 130                       + std::to_string(k) + std::string("] and x[")
 131                       + std::to_string(i) + std::string("] is ")
 132                       + std::string("smaller than the epsilon of ")
 133                       + std::string(typeid(Real).name());
 134                     throw std::logic_error(msg);
 135                 }
 136                 inv_product *= diff;
 137             }
 138             if (i % 2 == 0)
 139             {
 140                 m_w[k] += 1/inv_product;
 141             }
 142             else
 143             {
 144                 m_w[k] -= 1/inv_product;
 145             }
 146         }
 147     }
 148 }
 149
 150
 151 template<class Real>
 152 Real barycentric_rational_imp<Real>::operator()(Real x) const
 153 {
 154     Real numerator = 0;
 155     Real denominator = 0;
 156     for(size_t i = 0; i < m_x.size(); ++i)
 157     {
 158         // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
 159         // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
 160         // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
 161         if (x == m_x[i])
 162         {
 163             return m_y[i];
 164         }
 165         Real t = m_w[i]/(x - m_x[i]);
 166         numerator += t*m_y[i];
 167         denominator += t;
 168     }
 169     return numerator/denominator;
 170 }
 171
 172 /*
 173  * A formula for computing the derivative of the barycentric representation is given in
 174  * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
 175  * Mathematics of Computation, v47, number 175, 1986.
 176  * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
 177  * and reviewed in
 178  * Recent developments in barycentric rational interpolation
 179  * Jean-Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
 180  *
 181  * Is it possible to complete this in one pass through the data?
 182  */
 183
 184 template<class Real>
 185 Real barycentric_rational_imp<Real>::prime(Real x) const
 186 {
 187     Real rx = this->operator()(x);
 188     Real numerator = 0;
 189     Real denominator = 0;
 190     for(size_t i = 0; i < m_x.size(); ++i)
 191     {
 192         if (x == m_x[i])
 193         {
 194             Real sum = 0;
 195             for (size_t j = 0; j < m_x.size(); ++j)
 196             {
 197                 if (j == i)
 198                 {
 199                     continue;
 200                 }
 201                 sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]);
 202             }
 203             return -sum/m_w[i];
 204         }
 205         Real t = m_w[i]/(x - m_x[i]);
 206         Real diff = (rx - m_y[i])/(x-m_x[i]);
 207         numerator += t*diff;
 208         denominator += t;
 209     }
 210
 211     return numerator/denominator;
 212 }
 213 }}}}
 214 #endif