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