]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
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 <boost/lexical_cast.hpp> | |
13 | #include <boost/math/special_functions/fpclassify.hpp> | |
14 | #include <boost/core/demangle.hpp> | |
15 | ||
16 | namespace boost{ namespace math{ namespace detail{ | |
17 | ||
18 | template<class Real> | |
19 | class barycentric_rational_imp | |
20 | { | |
21 | public: | |
22 | template <class InputIterator1, class InputIterator2> | |
23 | barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3); | |
24 | ||
25 | Real operator()(Real x) const; | |
26 | ||
27 | // The barycentric weights are not really that interesting; except to the unit tests! | |
28 | Real weight(size_t i) const { return m_w[i]; } | |
29 | ||
30 | private: | |
31 | // Technically, we do not need to copy over y to m_y, or x to m_x. | |
32 | // We could simply store a pointer. However, in doing so, | |
33 | // we'd need to make sure the user managed the lifetime of m_y, | |
34 | // and didn't alter its data. Since we are unlikely to run out of | |
35 | // memory for a linearly scaling algorithm, it seems best just to make a copy. | |
36 | std::vector<Real> m_y; | |
37 | std::vector<Real> m_x; | |
38 | std::vector<Real> m_w; | |
39 | }; | |
40 | ||
41 | template <class Real> | |
42 | template <class InputIterator1, class InputIterator2> | |
43 | barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order) | |
44 | { | |
45 | using std::abs; | |
46 | ||
47 | std::ptrdiff_t n = std::distance(start_x, end_x); | |
48 | ||
49 | if (approximation_order >= (std::size_t)n) | |
50 | { | |
51 | throw std::domain_error("Approximation order must be < data length."); | |
52 | } | |
53 | ||
54 | // Big sad memcpy to make sure the object is easy to use. | |
55 | m_x.resize(n); | |
56 | m_y.resize(n); | |
57 | for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i) | |
58 | { | |
59 | // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy: | |
60 | if(boost::math::isnan(*start_x)) | |
61 | { | |
62 | std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN"; | |
63 | throw std::domain_error(msg); | |
64 | } | |
65 | ||
66 | if(boost::math::isnan(*start_y)) | |
67 | { | |
68 | std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN"; | |
69 | throw std::domain_error(msg); | |
70 | } | |
71 | ||
72 | m_x[i] = *start_x; | |
73 | m_y[i] = *start_y; | |
74 | } | |
75 | ||
76 | m_w.resize(n, 0); | |
77 | for(int64_t k = 0; k < n; ++k) | |
78 | { | |
79 | int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0); | |
80 | int64_t i_max = k; | |
81 | if (k >= n - (std::ptrdiff_t)approximation_order) | |
82 | { | |
83 | i_max = n - approximation_order - 1; | |
84 | } | |
85 | ||
86 | for(int64_t i = i_min; i <= i_max; ++i) | |
87 | { | |
88 | Real inv_product = 1; | |
89 | int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1)); | |
90 | for(int64_t j = i; j <= j_max; ++j) | |
91 | { | |
92 | if (j == k) | |
93 | { | |
94 | continue; | |
95 | } | |
96 | ||
97 | Real diff = m_x[k] - m_x[j]; | |
98 | if (abs(diff) < std::numeric_limits<Real>::epsilon()) | |
99 | { | |
100 | std::string msg = std::string("Spacing between x[") | |
101 | + boost::lexical_cast<std::string>(k) + std::string("] and x[") | |
102 | + boost::lexical_cast<std::string>(i) + std::string("] is ") | |
103 | + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ") | |
104 | + boost::core::demangle(typeid(Real).name()); | |
105 | throw std::logic_error(msg); | |
106 | } | |
107 | inv_product *= diff; | |
108 | } | |
109 | if (i % 2 == 0) | |
110 | { | |
111 | m_w[k] += 1/inv_product; | |
112 | } | |
113 | else | |
114 | { | |
115 | m_w[k] -= 1/inv_product; | |
116 | } | |
117 | } | |
118 | } | |
119 | } | |
120 | ||
121 | template<class Real> | |
122 | Real barycentric_rational_imp<Real>::operator()(Real x) const | |
123 | { | |
124 | Real numerator = 0; | |
125 | Real denominator = 0; | |
126 | for(size_t i = 0; i < m_x.size(); ++i) | |
127 | { | |
128 | // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality. | |
129 | // 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, | |
130 | // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715 | |
131 | if (x == m_x[i]) | |
132 | { | |
133 | return m_y[i]; | |
134 | } | |
135 | Real t = m_w[i]/(x - m_x[i]); | |
136 | numerator += t*m_y[i]; | |
137 | denominator += t; | |
138 | } | |
139 | return numerator/denominator; | |
140 | } | |
141 | ||
142 | /* | |
143 | * A formula for computing the derivative of the barycentric representation is given in | |
144 | * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner, | |
145 | * Mathematics of Computation, v47, number 175, 1986. | |
146 | * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf | |
147 | * However, this requires a lot of machinery which is not built into the library at present. | |
148 | * So we wait until there is a requirement to interpolate the derivative. | |
149 | */ | |
150 | }}} | |
151 | #endif |