| 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 | |