boost/math/tools/rational.hpp

   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);
 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*)
 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)
 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)
 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)
 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)
 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)
 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)
 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)
 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)
 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)
 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);
 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*)
 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)
 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)
 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)
 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