]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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_HEADER() <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_HEADER() <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 mpl::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 mpl::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 mpl::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 mpl::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 mpl::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<mpl::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 |