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