math/tools/rational.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_RATIONAL_HPP
0007 #define BOOST_MATH_TOOLS_RATIONAL_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/assert.hpp>
0015 #include <boost/math/tools/type_traits.hpp>
0016 #include <boost/math/tools/cstdint.hpp>
0017
0018 #ifndef BOOST_MATH_HAS_NVRTC
0019 #include <array>
0020 #endif
0021
0022 #if BOOST_MATH_POLY_METHOD == 1
0023 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/polynomial_horner1_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0024 #  include BOOST_HEADER()
0025 #  undef BOOST_HEADER
0026 #elif BOOST_MATH_POLY_METHOD == 2
0027 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/polynomial_horner2_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0028 #  include BOOST_HEADER()
0029 #  undef BOOST_HEADER
0030 #elif BOOST_MATH_POLY_METHOD == 3
0031 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/polynomial_horner3_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0032 #  include BOOST_HEADER()
0033 #  undef BOOST_HEADER
0034 #endif
0035 #if BOOST_MATH_RATIONAL_METHOD == 1
0036 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/rational_horner1_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0037 #  include BOOST_HEADER()
0038 #  undef BOOST_HEADER
0039 #elif BOOST_MATH_RATIONAL_METHOD == 2
0040 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/rational_horner2_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0041 #  include BOOST_HEADER()
0042 #  undef BOOST_HEADER
0043 #elif BOOST_MATH_RATIONAL_METHOD == 3
0044 #  define BOOST_HEADER() <BOOST_MATH_JOIN(boost/math/tools/detail/rational_horner3_, BOOST_MATH_MAX_POLY_ORDER).hpp>
0045 #  include BOOST_HEADER()
0046 #  undef BOOST_HEADER
0047 #endif
0048
0049 #if 0
0050 //
0051 // This just allows dependency trackers to find the headers
0052 // used in the above PP-magic.
0053 //
0054 #include <boost/math/tools/detail/polynomial_horner1_2.hpp>
0055 #include <boost/math/tools/detail/polynomial_horner1_3.hpp>
0056 #include <boost/math/tools/detail/polynomial_horner1_4.hpp>
0057 #include <boost/math/tools/detail/polynomial_horner1_5.hpp>
0058 #include <boost/math/tools/detail/polynomial_horner1_6.hpp>
0059 #include <boost/math/tools/detail/polynomial_horner1_7.hpp>
0060 #include <boost/math/tools/detail/polynomial_horner1_8.hpp>
0061 #include <boost/math/tools/detail/polynomial_horner1_9.hpp>
0062 #include <boost/math/tools/detail/polynomial_horner1_10.hpp>
0063 #include <boost/math/tools/detail/polynomial_horner1_11.hpp>
0064 #include <boost/math/tools/detail/polynomial_horner1_12.hpp>
0065 #include <boost/math/tools/detail/polynomial_horner1_13.hpp>
0066 #include <boost/math/tools/detail/polynomial_horner1_14.hpp>
0067 #include <boost/math/tools/detail/polynomial_horner1_15.hpp>
0068 #include <boost/math/tools/detail/polynomial_horner1_16.hpp>
0069 #include <boost/math/tools/detail/polynomial_horner1_17.hpp>
0070 #include <boost/math/tools/detail/polynomial_horner1_18.hpp>
0071 #include <boost/math/tools/detail/polynomial_horner1_19.hpp>
0072 #include <boost/math/tools/detail/polynomial_horner1_20.hpp>
0073 #include <boost/math/tools/detail/polynomial_horner2_2.hpp>
0074 #include <boost/math/tools/detail/polynomial_horner2_3.hpp>
0075 #include <boost/math/tools/detail/polynomial_horner2_4.hpp>
0076 #include <boost/math/tools/detail/polynomial_horner2_5.hpp>
0077 #include <boost/math/tools/detail/polynomial_horner2_6.hpp>
0078 #include <boost/math/tools/detail/polynomial_horner2_7.hpp>
0079 #include <boost/math/tools/detail/polynomial_horner2_8.hpp>
0080 #include <boost/math/tools/detail/polynomial_horner2_9.hpp>
0081 #include <boost/math/tools/detail/polynomial_horner2_10.hpp>
0082 #include <boost/math/tools/detail/polynomial_horner2_11.hpp>
0083 #include <boost/math/tools/detail/polynomial_horner2_12.hpp>
0084 #include <boost/math/tools/detail/polynomial_horner2_13.hpp>
0085 #include <boost/math/tools/detail/polynomial_horner2_14.hpp>
0086 #include <boost/math/tools/detail/polynomial_horner2_15.hpp>
0087 #include <boost/math/tools/detail/polynomial_horner2_16.hpp>
0088 #include <boost/math/tools/detail/polynomial_horner2_17.hpp>
0089 #include <boost/math/tools/detail/polynomial_horner2_18.hpp>
0090 #include <boost/math/tools/detail/polynomial_horner2_19.hpp>
0091 #include <boost/math/tools/detail/polynomial_horner2_20.hpp>
0092 #include <boost/math/tools/detail/polynomial_horner3_2.hpp>
0093 #include <boost/math/tools/detail/polynomial_horner3_3.hpp>
0094 #include <boost/math/tools/detail/polynomial_horner3_4.hpp>
0095 #include <boost/math/tools/detail/polynomial_horner3_5.hpp>
0096 #include <boost/math/tools/detail/polynomial_horner3_6.hpp>
0097 #include <boost/math/tools/detail/polynomial_horner3_7.hpp>
0098 #include <boost/math/tools/detail/polynomial_horner3_8.hpp>
0099 #include <boost/math/tools/detail/polynomial_horner3_9.hpp>
0100 #include <boost/math/tools/detail/polynomial_horner3_10.hpp>
0101 #include <boost/math/tools/detail/polynomial_horner3_11.hpp>
0102 #include <boost/math/tools/detail/polynomial_horner3_12.hpp>
0103 #include <boost/math/tools/detail/polynomial_horner3_13.hpp>
0104 #include <boost/math/tools/detail/polynomial_horner3_14.hpp>
0105 #include <boost/math/tools/detail/polynomial_horner3_15.hpp>
0106 #include <boost/math/tools/detail/polynomial_horner3_16.hpp>
0107 #include <boost/math/tools/detail/polynomial_horner3_17.hpp>
0108 #include <boost/math/tools/detail/polynomial_horner3_18.hpp>
0109 #include <boost/math/tools/detail/polynomial_horner3_19.hpp>
0110 #include <boost/math/tools/detail/polynomial_horner3_20.hpp>
0111 #include <boost/math/tools/detail/rational_horner1_2.hpp>
0112 #include <boost/math/tools/detail/rational_horner1_3.hpp>
0113 #include <boost/math/tools/detail/rational_horner1_4.hpp>
0114 #include <boost/math/tools/detail/rational_horner1_5.hpp>
0115 #include <boost/math/tools/detail/rational_horner1_6.hpp>
0116 #include <boost/math/tools/detail/rational_horner1_7.hpp>
0117 #include <boost/math/tools/detail/rational_horner1_8.hpp>
0118 #include <boost/math/tools/detail/rational_horner1_9.hpp>
0119 #include <boost/math/tools/detail/rational_horner1_10.hpp>
0120 #include <boost/math/tools/detail/rational_horner1_11.hpp>
0121 #include <boost/math/tools/detail/rational_horner1_12.hpp>
0122 #include <boost/math/tools/detail/rational_horner1_13.hpp>
0123 #include <boost/math/tools/detail/rational_horner1_14.hpp>
0124 #include <boost/math/tools/detail/rational_horner1_15.hpp>
0125 #include <boost/math/tools/detail/rational_horner1_16.hpp>
0126 #include <boost/math/tools/detail/rational_horner1_17.hpp>
0127 #include <boost/math/tools/detail/rational_horner1_18.hpp>
0128 #include <boost/math/tools/detail/rational_horner1_19.hpp>
0129 #include <boost/math/tools/detail/rational_horner1_20.hpp>
0130 #include <boost/math/tools/detail/rational_horner2_2.hpp>
0131 #include <boost/math/tools/detail/rational_horner2_3.hpp>
0132 #include <boost/math/tools/detail/rational_horner2_4.hpp>
0133 #include <boost/math/tools/detail/rational_horner2_5.hpp>
0134 #include <boost/math/tools/detail/rational_horner2_6.hpp>
0135 #include <boost/math/tools/detail/rational_horner2_7.hpp>
0136 #include <boost/math/tools/detail/rational_horner2_8.hpp>
0137 #include <boost/math/tools/detail/rational_horner2_9.hpp>
0138 #include <boost/math/tools/detail/rational_horner2_10.hpp>
0139 #include <boost/math/tools/detail/rational_horner2_11.hpp>
0140 #include <boost/math/tools/detail/rational_horner2_12.hpp>
0141 #include <boost/math/tools/detail/rational_horner2_13.hpp>
0142 #include <boost/math/tools/detail/rational_horner2_14.hpp>
0143 #include <boost/math/tools/detail/rational_horner2_15.hpp>
0144 #include <boost/math/tools/detail/rational_horner2_16.hpp>
0145 #include <boost/math/tools/detail/rational_horner2_17.hpp>
0146 #include <boost/math/tools/detail/rational_horner2_18.hpp>
0147 #include <boost/math/tools/detail/rational_horner2_19.hpp>
0148 #include <boost/math/tools/detail/rational_horner2_20.hpp>
0149 #include <boost/math/tools/detail/rational_horner3_2.hpp>
0150 #include <boost/math/tools/detail/rational_horner3_3.hpp>
0151 #include <boost/math/tools/detail/rational_horner3_4.hpp>
0152 #include <boost/math/tools/detail/rational_horner3_5.hpp>
0153 #include <boost/math/tools/detail/rational_horner3_6.hpp>
0154 #include <boost/math/tools/detail/rational_horner3_7.hpp>
0155 #include <boost/math/tools/detail/rational_horner3_8.hpp>
0156 #include <boost/math/tools/detail/rational_horner3_9.hpp>
0157 #include <boost/math/tools/detail/rational_horner3_10.hpp>
0158 #include <boost/math/tools/detail/rational_horner3_11.hpp>
0159 #include <boost/math/tools/detail/rational_horner3_12.hpp>
0160 #include <boost/math/tools/detail/rational_horner3_13.hpp>
0161 #include <boost/math/tools/detail/rational_horner3_14.hpp>
0162 #include <boost/math/tools/detail/rational_horner3_15.hpp>
0163 #include <boost/math/tools/detail/rational_horner3_16.hpp>
0164 #include <boost/math/tools/detail/rational_horner3_17.hpp>
0165 #include <boost/math/tools/detail/rational_horner3_18.hpp>
0166 #include <boost/math/tools/detail/rational_horner3_19.hpp>
0167 #include <boost/math/tools/detail/rational_horner3_20.hpp>
0168 #endif
0169
0170 namespace boost{ namespace math{ namespace tools{
0171
0172 //
0173 // Forward declaration to keep two phase lookup happy:
0174 //
0175 template <class T, class U>
0176 BOOST_MATH_GPU_ENABLED U evaluate_polynomial(const T* poly, U const& z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U);
0177
0178 namespace detail{
0179
0180 template <class T, class V, class Tag>
0181 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_polynomial_c_imp(const T* a, const V& val, const Tag*) BOOST_MATH_NOEXCEPT(V)
0182 {
0183    return evaluate_polynomial(a, val, Tag::value);
0184 }
0185
0186 } // namespace detail
0187
0188 //
0189 // Polynomial evaluation with runtime size.
0190 // This requires a for-loop which may be more expensive than
0191 // the loop expanded versions above:
0192 //
0193 template <class T, class U>
0194 BOOST_MATH_GPU_ENABLED inline U evaluate_polynomial(const T* poly, U const& z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
0195 {
0196    BOOST_MATH_ASSERT(count > 0);
0197    U sum = static_cast<U>(poly[count - 1]);
0198    for(int i = static_cast<int>(count) - 2; i >= 0; --i)
0199    {
0200       sum *= z;
0201       sum += static_cast<U>(poly[i]);
0202    }
0203    return sum;
0204 }
0205 //
0206 // Compile time sized polynomials, just inline forwarders to the
0207 // implementations above:
0208 //
0209 template <boost::math::size_t N, class T, class V>
0210 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_polynomial(const T(&a)[N], const V& val) BOOST_MATH_NOEXCEPT(V)
0211 {
0212    typedef boost::math::integral_constant<int, N> tag_type;
0213    return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(nullptr));
0214 }
0215
0216 #ifndef BOOST_MATH_HAS_NVRTC
0217 template <boost::math::size_t N, class T, class V>
0218 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_polynomial(const std::array<T,N>& a, const V& val) BOOST_MATH_NOEXCEPT(V)
0219 {
0220    typedef boost::math::integral_constant<int, N> tag_type;
0221    return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()), val, static_cast<tag_type const*>(nullptr));
0222 }
0223 #endif
0224 //
0225 // Even polynomials are trivial: just square the argument!
0226 //
0227 template <class T, class U>
0228 BOOST_MATH_GPU_ENABLED inline U evaluate_even_polynomial(const T* poly, U z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
0229 {
0230    return evaluate_polynomial(poly, U(z*z), count);
0231 }
0232
0233 template <boost::math::size_t N, class T, class V>
0234 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_even_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V)
0235 {
0236    return evaluate_polynomial(a, V(z*z));
0237 }
0238
0239 #ifndef BOOST_MATH_HAS_NVRTC
0240 template <boost::math::size_t N, class T, class V>
0241 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_even_polynomial(const std::array<T,N>& a, const V& z) BOOST_MATH_NOEXCEPT(V)
0242 {
0243    return evaluate_polynomial(a, V(z*z));
0244 }
0245 #endif
0246 //
0247 // Odd polynomials come next:
0248 //
0249 template <class T, class U>
0250 BOOST_MATH_GPU_ENABLED inline U evaluate_odd_polynomial(const T* poly, U z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
0251 {
0252    return poly[0] + z * evaluate_polynomial(poly+1, U(z*z), count-1);
0253 }
0254
0255 template <boost::math::size_t N, class T, class V>
0256 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_odd_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V)
0257 {
0258    typedef boost::math::integral_constant<int, N-1> tag_type;
0259    return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, V(z*z), static_cast<tag_type const*>(nullptr));
0260 }
0261
0262 #ifndef BOOST_MATH_HAS_NVRTC
0263 template <boost::math::size_t N, class T, class V>
0264 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_odd_polynomial(const std::array<T,N>& a, const V& z) BOOST_MATH_NOEXCEPT(V)
0265 {
0266    typedef boost::math::integral_constant<int, N-1> tag_type;
0267    return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, V(z*z), static_cast<tag_type const*>(nullptr));
0268 }
0269 #endif
0270
0271 template <class T, class U, class V>
0272 BOOST_MATH_GPU_ENABLED V evaluate_rational(const T* num, const U* denom, const V& z_, boost::math::size_t count) BOOST_MATH_NOEXCEPT(V);
0273
0274 namespace detail{
0275
0276 template <class T, class U, class V, class Tag>
0277 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_rational_c_imp(const T* num, const U* denom, const V& z, const Tag*) BOOST_MATH_NOEXCEPT(V)
0278 {
0279    return boost::math::tools::evaluate_rational(num, denom, z, Tag::value);
0280 }
0281
0282 }
0283 //
0284 // Rational functions: numerator and denominator must be
0285 // equal in size.  These always have a for-loop and so may be less
0286 // efficient than evaluating a pair of polynomials. However, there
0287 // are some tricks we can use to prevent overflow that might otherwise
0288 // occur in polynomial evaluation, if z is large.  This is important
0289 // in our Lanczos code for example.
0290 //
0291 template <class T, class U, class V>
0292 BOOST_MATH_GPU_ENABLED V evaluate_rational(const T* num, const U* denom, const V& z_, boost::math::size_t count) BOOST_MATH_NOEXCEPT(V)
0293 {
0294    V z(z_);
0295    V s1, s2;
0296    if(z <= 1)
0297    {
0298       s1 = static_cast<V>(num[count-1]);
0299       s2 = static_cast<V>(denom[count-1]);
0300       for(int i = (int)count - 2; i >= 0; --i)
0301       {
0302          s1 *= z;
0303          s2 *= z;
0304          s1 += num[i];
0305          s2 += denom[i];
0306       }
0307    }
0308    else
0309    {
0310       z = 1 / z;
0311       s1 = static_cast<V>(num[0]);
0312       s2 = static_cast<V>(denom[0]);
0313       for(unsigned i = 1; i < count; ++i)
0314       {
0315          s1 *= z;
0316          s2 *= z;
0317          s1 += num[i];
0318          s2 += denom[i];
0319       }
0320    }
0321    return s1 / s2;
0322 }
0323
0324 template <boost::math::size_t N, class T, class U, class V>
0325 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_rational(const T(&a)[N], const U(&b)[N], const V& z) BOOST_MATH_NOEXCEPT(V)
0326 {
0327    return detail::evaluate_rational_c_imp(a, b, z, static_cast<const boost::math::integral_constant<int, N>*>(nullptr));
0328 }
0329
0330 #ifndef BOOST_MATH_HAS_NVRTC
0331 template <boost::math::size_t N, class T, class U, class V>
0332 BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_rational(const std::array<T,N>& a, const std::array<U,N>& b, const V& z) BOOST_MATH_NOEXCEPT(V)
0333 {
0334    return detail::evaluate_rational_c_imp(a.data(), b.data(), z, static_cast<boost::math::integral_constant<int, N>*>(nullptr));
0335 }
0336 #endif
0337
0338 } // namespace tools
0339 } // namespace math
0340 } // namespace boost
0341
0342 #endif // BOOST_MATH_TOOLS_RATIONAL_HPP
0343
0344
0345
0346