1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the |
3 | // Boost Software License, Version 1.0. (See accompanying file |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
5 | |
6 | #ifndef BOOST_MATH_TOOLS_RATIONAL_HPP |
7 | #define BOOST_MATH_TOOLS_RATIONAL_HPP |
8 | |
9 | #ifdef _MSC_VER |
10 | #pragma once |
11 | #endif |
12 | |
13 | #include <boost/array.hpp> |
14 | #include <boost/math/tools/config.hpp> |
15 | #include <boost/mpl/int.hpp> |
16 | |
17 | #if BOOST_MATH_POLY_METHOD == 1 |
18 | # define BOOST_HEADER() <BOOST_JOIN(boost/math/tools/detail/polynomial_horner1_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
19 | # include BOOST_HEADER() |
20 | # undef BOOST_HEADER |
21 | #elif BOOST_MATH_POLY_METHOD == 2 |
22 | # define BOOST_HEADER() <BOOST_JOIN(boost/math/tools/detail/polynomial_horner2_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
23 | # include BOOST_HEADER() |
24 | # undef BOOST_HEADER |
25 | #elif BOOST_MATH_POLY_METHOD == 3 |
26 | # define () <BOOST_JOIN(boost/math/tools/detail/polynomial_horner3_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
27 | # include BOOST_HEADER() |
28 | # undef BOOST_HEADER |
29 | #endif |
30 | #if BOOST_MATH_RATIONAL_METHOD == 1 |
31 | # define BOOST_HEADER() <BOOST_JOIN(boost/math/tools/detail/rational_horner1_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
32 | # include BOOST_HEADER() |
33 | # undef BOOST_HEADER |
34 | #elif BOOST_MATH_RATIONAL_METHOD == 2 |
35 | # define BOOST_HEADER() <BOOST_JOIN(boost/math/tools/detail/rational_horner2_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
36 | # include BOOST_HEADER() |
37 | # undef BOOST_HEADER |
38 | #elif BOOST_MATH_RATIONAL_METHOD == 3 |
39 | # define () <BOOST_JOIN(boost/math/tools/detail/rational_horner3_, BOOST_MATH_MAX_POLY_ORDER).hpp> |
40 | # include BOOST_HEADER() |
41 | # undef BOOST_HEADER |
42 | #endif |
43 | |
44 | #if 0 |
45 | // |
46 | // This just allows dependency trackers to find the headers |
47 | // used in the above PP-magic. |
48 | // |
49 | #include <boost/math/tools/detail/polynomial_horner1_2.hpp> |
50 | #include <boost/math/tools/detail/polynomial_horner1_3.hpp> |
51 | #include <boost/math/tools/detail/polynomial_horner1_4.hpp> |
52 | #include <boost/math/tools/detail/polynomial_horner1_5.hpp> |
53 | #include <boost/math/tools/detail/polynomial_horner1_6.hpp> |
54 | #include <boost/math/tools/detail/polynomial_horner1_7.hpp> |
55 | #include <boost/math/tools/detail/polynomial_horner1_8.hpp> |
56 | #include <boost/math/tools/detail/polynomial_horner1_9.hpp> |
57 | #include <boost/math/tools/detail/polynomial_horner1_10.hpp> |
58 | #include <boost/math/tools/detail/polynomial_horner1_11.hpp> |
59 | #include <boost/math/tools/detail/polynomial_horner1_12.hpp> |
60 | #include <boost/math/tools/detail/polynomial_horner1_13.hpp> |
61 | #include <boost/math/tools/detail/polynomial_horner1_14.hpp> |
62 | #include <boost/math/tools/detail/polynomial_horner1_15.hpp> |
63 | #include <boost/math/tools/detail/polynomial_horner1_16.hpp> |
64 | #include <boost/math/tools/detail/polynomial_horner1_17.hpp> |
65 | #include <boost/math/tools/detail/polynomial_horner1_18.hpp> |
66 | #include <boost/math/tools/detail/polynomial_horner1_19.hpp> |
67 | #include <boost/math/tools/detail/polynomial_horner1_20.hpp> |
68 | #include <boost/math/tools/detail/polynomial_horner2_2.hpp> |
69 | #include <boost/math/tools/detail/polynomial_horner2_3.hpp> |
70 | #include <boost/math/tools/detail/polynomial_horner2_4.hpp> |
71 | #include <boost/math/tools/detail/polynomial_horner2_5.hpp> |
72 | #include <boost/math/tools/detail/polynomial_horner2_6.hpp> |
73 | #include <boost/math/tools/detail/polynomial_horner2_7.hpp> |
74 | #include <boost/math/tools/detail/polynomial_horner2_8.hpp> |
75 | #include <boost/math/tools/detail/polynomial_horner2_9.hpp> |
76 | #include <boost/math/tools/detail/polynomial_horner2_10.hpp> |
77 | #include <boost/math/tools/detail/polynomial_horner2_11.hpp> |
78 | #include <boost/math/tools/detail/polynomial_horner2_12.hpp> |
79 | #include <boost/math/tools/detail/polynomial_horner2_13.hpp> |
80 | #include <boost/math/tools/detail/polynomial_horner2_14.hpp> |
81 | #include <boost/math/tools/detail/polynomial_horner2_15.hpp> |
82 | #include <boost/math/tools/detail/polynomial_horner2_16.hpp> |
83 | #include <boost/math/tools/detail/polynomial_horner2_17.hpp> |
84 | #include <boost/math/tools/detail/polynomial_horner2_18.hpp> |
85 | #include <boost/math/tools/detail/polynomial_horner2_19.hpp> |
86 | #include <boost/math/tools/detail/polynomial_horner2_20.hpp> |
87 | #include <boost/math/tools/detail/polynomial_horner3_2.hpp> |
88 | #include <boost/math/tools/detail/polynomial_horner3_3.hpp> |
89 | #include <boost/math/tools/detail/polynomial_horner3_4.hpp> |
90 | #include <boost/math/tools/detail/polynomial_horner3_5.hpp> |
91 | #include <boost/math/tools/detail/polynomial_horner3_6.hpp> |
92 | #include <boost/math/tools/detail/polynomial_horner3_7.hpp> |
93 | #include <boost/math/tools/detail/polynomial_horner3_8.hpp> |
94 | #include <boost/math/tools/detail/polynomial_horner3_9.hpp> |
95 | #include <boost/math/tools/detail/polynomial_horner3_10.hpp> |
96 | #include <boost/math/tools/detail/polynomial_horner3_11.hpp> |
97 | #include <boost/math/tools/detail/polynomial_horner3_12.hpp> |
98 | #include <boost/math/tools/detail/polynomial_horner3_13.hpp> |
99 | #include <boost/math/tools/detail/polynomial_horner3_14.hpp> |
100 | #include <boost/math/tools/detail/polynomial_horner3_15.hpp> |
101 | #include <boost/math/tools/detail/polynomial_horner3_16.hpp> |
102 | #include <boost/math/tools/detail/polynomial_horner3_17.hpp> |
103 | #include <boost/math/tools/detail/polynomial_horner3_18.hpp> |
104 | #include <boost/math/tools/detail/polynomial_horner3_19.hpp> |
105 | #include <boost/math/tools/detail/polynomial_horner3_20.hpp> |
106 | #include <boost/math/tools/detail/rational_horner1_2.hpp> |
107 | #include <boost/math/tools/detail/rational_horner1_3.hpp> |
108 | #include <boost/math/tools/detail/rational_horner1_4.hpp> |
109 | #include <boost/math/tools/detail/rational_horner1_5.hpp> |
110 | #include <boost/math/tools/detail/rational_horner1_6.hpp> |
111 | #include <boost/math/tools/detail/rational_horner1_7.hpp> |
112 | #include <boost/math/tools/detail/rational_horner1_8.hpp> |
113 | #include <boost/math/tools/detail/rational_horner1_9.hpp> |
114 | #include <boost/math/tools/detail/rational_horner1_10.hpp> |
115 | #include <boost/math/tools/detail/rational_horner1_11.hpp> |
116 | #include <boost/math/tools/detail/rational_horner1_12.hpp> |
117 | #include <boost/math/tools/detail/rational_horner1_13.hpp> |
118 | #include <boost/math/tools/detail/rational_horner1_14.hpp> |
119 | #include <boost/math/tools/detail/rational_horner1_15.hpp> |
120 | #include <boost/math/tools/detail/rational_horner1_16.hpp> |
121 | #include <boost/math/tools/detail/rational_horner1_17.hpp> |
122 | #include <boost/math/tools/detail/rational_horner1_18.hpp> |
123 | #include <boost/math/tools/detail/rational_horner1_19.hpp> |
124 | #include <boost/math/tools/detail/rational_horner1_20.hpp> |
125 | #include <boost/math/tools/detail/rational_horner2_2.hpp> |
126 | #include <boost/math/tools/detail/rational_horner2_3.hpp> |
127 | #include <boost/math/tools/detail/rational_horner2_4.hpp> |
128 | #include <boost/math/tools/detail/rational_horner2_5.hpp> |
129 | #include <boost/math/tools/detail/rational_horner2_6.hpp> |
130 | #include <boost/math/tools/detail/rational_horner2_7.hpp> |
131 | #include <boost/math/tools/detail/rational_horner2_8.hpp> |
132 | #include <boost/math/tools/detail/rational_horner2_9.hpp> |
133 | #include <boost/math/tools/detail/rational_horner2_10.hpp> |
134 | #include <boost/math/tools/detail/rational_horner2_11.hpp> |
135 | #include <boost/math/tools/detail/rational_horner2_12.hpp> |
136 | #include <boost/math/tools/detail/rational_horner2_13.hpp> |
137 | #include <boost/math/tools/detail/rational_horner2_14.hpp> |
138 | #include <boost/math/tools/detail/rational_horner2_15.hpp> |
139 | #include <boost/math/tools/detail/rational_horner2_16.hpp> |
140 | #include <boost/math/tools/detail/rational_horner2_17.hpp> |
141 | #include <boost/math/tools/detail/rational_horner2_18.hpp> |
142 | #include <boost/math/tools/detail/rational_horner2_19.hpp> |
143 | #include <boost/math/tools/detail/rational_horner2_20.hpp> |
144 | #include <boost/math/tools/detail/rational_horner3_2.hpp> |
145 | #include <boost/math/tools/detail/rational_horner3_3.hpp> |
146 | #include <boost/math/tools/detail/rational_horner3_4.hpp> |
147 | #include <boost/math/tools/detail/rational_horner3_5.hpp> |
148 | #include <boost/math/tools/detail/rational_horner3_6.hpp> |
149 | #include <boost/math/tools/detail/rational_horner3_7.hpp> |
150 | #include <boost/math/tools/detail/rational_horner3_8.hpp> |
151 | #include <boost/math/tools/detail/rational_horner3_9.hpp> |
152 | #include <boost/math/tools/detail/rational_horner3_10.hpp> |
153 | #include <boost/math/tools/detail/rational_horner3_11.hpp> |
154 | #include <boost/math/tools/detail/rational_horner3_12.hpp> |
155 | #include <boost/math/tools/detail/rational_horner3_13.hpp> |
156 | #include <boost/math/tools/detail/rational_horner3_14.hpp> |
157 | #include <boost/math/tools/detail/rational_horner3_15.hpp> |
158 | #include <boost/math/tools/detail/rational_horner3_16.hpp> |
159 | #include <boost/math/tools/detail/rational_horner3_17.hpp> |
160 | #include <boost/math/tools/detail/rational_horner3_18.hpp> |
161 | #include <boost/math/tools/detail/rational_horner3_19.hpp> |
162 | #include <boost/math/tools/detail/rational_horner3_20.hpp> |
163 | #endif |
164 | |
165 | namespace boost{ namespace math{ namespace tools{ |
166 | |
167 | // |
168 | // Forward declaration to keep two phase lookup happy: |
169 | // |
170 | template <class T, class U> |
171 | U evaluate_polynomial(const T* poly, U const& z, std::size_t count) BOOST_MATH_NOEXCEPT(U); |
172 | |
173 | namespace detail{ |
174 | |
175 | template <class T, class V, class Tag> |
176 | inline V evaluate_polynomial_c_imp(const T* a, const V& val, const Tag*) BOOST_MATH_NOEXCEPT(V) |
177 | { |
178 | return evaluate_polynomial(a, val, Tag::value); |
179 | } |
180 | |
181 | } // namespace detail |
182 | |
183 | // |
184 | // Polynomial evaluation with runtime size. |
185 | // This requires a for-loop which may be more expensive than |
186 | // the loop expanded versions above: |
187 | // |
188 | template <class T, class U> |
189 | inline U evaluate_polynomial(const T* poly, U const& z, std::size_t count) BOOST_MATH_NOEXCEPT(U) |
190 | { |
191 | BOOST_ASSERT(count > 0); |
192 | U sum = static_cast<U>(poly[count - 1]); |
193 | for(int i = static_cast<int>(count) - 2; i >= 0; --i) |
194 | { |
195 | sum *= z; |
196 | sum += static_cast<U>(poly[i]); |
197 | } |
198 | return sum; |
199 | } |
200 | // |
201 | // Compile time sized polynomials, just inline forwarders to the |
202 | // implementations above: |
203 | // |
204 | template <std::size_t N, class T, class V> |
205 | inline V evaluate_polynomial(const T(&a)[N], const V& val) BOOST_MATH_NOEXCEPT(V) |
206 | { |
207 | typedef boost::integral_constant<int, N> tag_type; |
208 | return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(0)); |
209 | } |
210 | |
211 | template <std::size_t N, class T, class V> |
212 | inline V evaluate_polynomial(const boost::array<T,N>& a, const V& val) BOOST_MATH_NOEXCEPT(V) |
213 | { |
214 | typedef boost::integral_constant<int, N> tag_type; |
215 | return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()), val, static_cast<tag_type const*>(0)); |
216 | } |
217 | // |
218 | // Even polynomials are trivial: just square the argument! |
219 | // |
220 | template <class T, class U> |
221 | inline U evaluate_even_polynomial(const T* poly, U z, std::size_t count) BOOST_MATH_NOEXCEPT(U) |
222 | { |
223 | return evaluate_polynomial(poly, U(z*z), count); |
224 | } |
225 | |
226 | template <std::size_t N, class T, class V> |
227 | inline V evaluate_even_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V) |
228 | { |
229 | return evaluate_polynomial(a, V(z*z)); |
230 | } |
231 | |
232 | template <std::size_t N, class T, class V> |
233 | inline V evaluate_even_polynomial(const boost::array<T,N>& a, const V& z) BOOST_MATH_NOEXCEPT(V) |
234 | { |
235 | return evaluate_polynomial(a, V(z*z)); |
236 | } |
237 | // |
238 | // Odd polynomials come next: |
239 | // |
240 | template <class T, class U> |
241 | inline U evaluate_odd_polynomial(const T* poly, U z, std::size_t count) BOOST_MATH_NOEXCEPT(U) |
242 | { |
243 | return poly[0] + z * evaluate_polynomial(poly+1, U(z*z), count-1); |
244 | } |
245 | |
246 | template <std::size_t N, class T, class V> |
247 | inline V evaluate_odd_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V) |
248 | { |
249 | typedef boost::integral_constant<int, N-1> tag_type; |
250 | return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, V(z*z), static_cast<tag_type const*>(0)); |
251 | } |
252 | |
253 | template <std::size_t N, class T, class V> |
254 | inline V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z) BOOST_MATH_NOEXCEPT(V) |
255 | { |
256 | typedef boost::integral_constant<int, N-1> tag_type; |
257 | return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, V(z*z), static_cast<tag_type const*>(0)); |
258 | } |
259 | |
260 | template <class T, class U, class V> |
261 | V evaluate_rational(const T* num, const U* denom, const V& z_, std::size_t count) BOOST_MATH_NOEXCEPT(V); |
262 | |
263 | namespace detail{ |
264 | |
265 | template <class T, class U, class V, class Tag> |
266 | inline V evaluate_rational_c_imp(const T* num, const U* denom, const V& z, const Tag*) BOOST_MATH_NOEXCEPT(V) |
267 | { |
268 | return boost::math::tools::evaluate_rational(num, denom, z, Tag::value); |
269 | } |
270 | |
271 | } |
272 | // |
273 | // Rational functions: numerator and denominator must be |
274 | // equal in size. These always have a for-loop and so may be less |
275 | // efficient than evaluating a pair of polynomials. However, there |
276 | // are some tricks we can use to prevent overflow that might otherwise |
277 | // occur in polynomial evaluation, if z is large. This is important |
278 | // in our Lanczos code for example. |
279 | // |
280 | template <class T, class U, class V> |
281 | V evaluate_rational(const T* num, const U* denom, const V& z_, std::size_t count) BOOST_MATH_NOEXCEPT(V) |
282 | { |
283 | V z(z_); |
284 | V s1, s2; |
285 | if(z <= 1) |
286 | { |
287 | s1 = static_cast<V>(num[count-1]); |
288 | s2 = static_cast<V>(denom[count-1]); |
289 | for(int i = (int)count - 2; i >= 0; --i) |
290 | { |
291 | s1 *= z; |
292 | s2 *= z; |
293 | s1 += num[i]; |
294 | s2 += denom[i]; |
295 | } |
296 | } |
297 | else |
298 | { |
299 | z = 1 / z; |
300 | s1 = static_cast<V>(num[0]); |
301 | s2 = static_cast<V>(denom[0]); |
302 | for(unsigned i = 1; i < count; ++i) |
303 | { |
304 | s1 *= z; |
305 | s2 *= z; |
306 | s1 += num[i]; |
307 | s2 += denom[i]; |
308 | } |
309 | } |
310 | return s1 / s2; |
311 | } |
312 | |
313 | template <std::size_t N, class T, class U, class V> |
314 | inline V evaluate_rational(const T(&a)[N], const U(&b)[N], const V& z) BOOST_MATH_NOEXCEPT(V) |
315 | { |
316 | return detail::evaluate_rational_c_imp(a, b, z, static_cast<const boost::integral_constant<int, N>*>(0)); |
317 | } |
318 | |
319 | template <std::size_t N, class T, class U, class V> |
320 | inline V evaluate_rational(const boost::array<T,N>& a, const boost::array<U,N>& b, const V& z) BOOST_MATH_NOEXCEPT(V) |
321 | { |
322 | return detail::evaluate_rational_c_imp(a.data(), b.data(), z, static_cast<boost::integral_constant<int, N>*>(0)); |
323 | } |
324 | |
325 | } // namespace tools |
326 | } // namespace math |
327 | } // namespace boost |
328 | |
329 | #endif // BOOST_MATH_TOOLS_RATIONAL_HPP |
330 | |
331 | |
332 | |
333 | |
334 | |