rational.hpp 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347
  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. #ifndef BOOST_MATH_TOOLS_RATIONAL_HPP
  6. #define BOOST_MATH_TOOLS_RATIONAL_HPP
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/math/tools/config.hpp>
  11. #include <boost/math/tools/assert.hpp>
  12. #include <boost/math/tools/type_traits.hpp>
  13. #include <boost/math/tools/cstdint.hpp>
  14. #ifndef BOOST_MATH_HAS_NVRTC
  15. #include <array>
  16. #endif
  17. #if BOOST_MATH_POLY_METHOD == 1
  18. # define BOOST_HEADER() <BOOST_MATH_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_MATH_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_MATH_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_MATH_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_MATH_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_MATH_JOIN(boost/math/tools/detail/rational_horner3_, BOOST_MATH_MAX_POLY_ORDER).hpp>
  40. # include BOOST_HEADER()
  41. # undef BOOST_HEADER
  42. #endif
  43. #if 0
  44. //
  45. // This just allows dependency trackers to find the headers
  46. // used in the above PP-magic.
  47. //
  48. #include <boost/math/tools/detail/polynomial_horner1_2.hpp>
  49. #include <boost/math/tools/detail/polynomial_horner1_3.hpp>
  50. #include <boost/math/tools/detail/polynomial_horner1_4.hpp>
  51. #include <boost/math/tools/detail/polynomial_horner1_5.hpp>
  52. #include <boost/math/tools/detail/polynomial_horner1_6.hpp>
  53. #include <boost/math/tools/detail/polynomial_horner1_7.hpp>
  54. #include <boost/math/tools/detail/polynomial_horner1_8.hpp>
  55. #include <boost/math/tools/detail/polynomial_horner1_9.hpp>
  56. #include <boost/math/tools/detail/polynomial_horner1_10.hpp>
  57. #include <boost/math/tools/detail/polynomial_horner1_11.hpp>
  58. #include <boost/math/tools/detail/polynomial_horner1_12.hpp>
  59. #include <boost/math/tools/detail/polynomial_horner1_13.hpp>
  60. #include <boost/math/tools/detail/polynomial_horner1_14.hpp>
  61. #include <boost/math/tools/detail/polynomial_horner1_15.hpp>
  62. #include <boost/math/tools/detail/polynomial_horner1_16.hpp>
  63. #include <boost/math/tools/detail/polynomial_horner1_17.hpp>
  64. #include <boost/math/tools/detail/polynomial_horner1_18.hpp>
  65. #include <boost/math/tools/detail/polynomial_horner1_19.hpp>
  66. #include <boost/math/tools/detail/polynomial_horner1_20.hpp>
  67. #include <boost/math/tools/detail/polynomial_horner2_2.hpp>
  68. #include <boost/math/tools/detail/polynomial_horner2_3.hpp>
  69. #include <boost/math/tools/detail/polynomial_horner2_4.hpp>
  70. #include <boost/math/tools/detail/polynomial_horner2_5.hpp>
  71. #include <boost/math/tools/detail/polynomial_horner2_6.hpp>
  72. #include <boost/math/tools/detail/polynomial_horner2_7.hpp>
  73. #include <boost/math/tools/detail/polynomial_horner2_8.hpp>
  74. #include <boost/math/tools/detail/polynomial_horner2_9.hpp>
  75. #include <boost/math/tools/detail/polynomial_horner2_10.hpp>
  76. #include <boost/math/tools/detail/polynomial_horner2_11.hpp>
  77. #include <boost/math/tools/detail/polynomial_horner2_12.hpp>
  78. #include <boost/math/tools/detail/polynomial_horner2_13.hpp>
  79. #include <boost/math/tools/detail/polynomial_horner2_14.hpp>
  80. #include <boost/math/tools/detail/polynomial_horner2_15.hpp>
  81. #include <boost/math/tools/detail/polynomial_horner2_16.hpp>
  82. #include <boost/math/tools/detail/polynomial_horner2_17.hpp>
  83. #include <boost/math/tools/detail/polynomial_horner2_18.hpp>
  84. #include <boost/math/tools/detail/polynomial_horner2_19.hpp>
  85. #include <boost/math/tools/detail/polynomial_horner2_20.hpp>
  86. #include <boost/math/tools/detail/polynomial_horner3_2.hpp>
  87. #include <boost/math/tools/detail/polynomial_horner3_3.hpp>
  88. #include <boost/math/tools/detail/polynomial_horner3_4.hpp>
  89. #include <boost/math/tools/detail/polynomial_horner3_5.hpp>
  90. #include <boost/math/tools/detail/polynomial_horner3_6.hpp>
  91. #include <boost/math/tools/detail/polynomial_horner3_7.hpp>
  92. #include <boost/math/tools/detail/polynomial_horner3_8.hpp>
  93. #include <boost/math/tools/detail/polynomial_horner3_9.hpp>
  94. #include <boost/math/tools/detail/polynomial_horner3_10.hpp>
  95. #include <boost/math/tools/detail/polynomial_horner3_11.hpp>
  96. #include <boost/math/tools/detail/polynomial_horner3_12.hpp>
  97. #include <boost/math/tools/detail/polynomial_horner3_13.hpp>
  98. #include <boost/math/tools/detail/polynomial_horner3_14.hpp>
  99. #include <boost/math/tools/detail/polynomial_horner3_15.hpp>
  100. #include <boost/math/tools/detail/polynomial_horner3_16.hpp>
  101. #include <boost/math/tools/detail/polynomial_horner3_17.hpp>
  102. #include <boost/math/tools/detail/polynomial_horner3_18.hpp>
  103. #include <boost/math/tools/detail/polynomial_horner3_19.hpp>
  104. #include <boost/math/tools/detail/polynomial_horner3_20.hpp>
  105. #include <boost/math/tools/detail/rational_horner1_2.hpp>
  106. #include <boost/math/tools/detail/rational_horner1_3.hpp>
  107. #include <boost/math/tools/detail/rational_horner1_4.hpp>
  108. #include <boost/math/tools/detail/rational_horner1_5.hpp>
  109. #include <boost/math/tools/detail/rational_horner1_6.hpp>
  110. #include <boost/math/tools/detail/rational_horner1_7.hpp>
  111. #include <boost/math/tools/detail/rational_horner1_8.hpp>
  112. #include <boost/math/tools/detail/rational_horner1_9.hpp>
  113. #include <boost/math/tools/detail/rational_horner1_10.hpp>
  114. #include <boost/math/tools/detail/rational_horner1_11.hpp>
  115. #include <boost/math/tools/detail/rational_horner1_12.hpp>
  116. #include <boost/math/tools/detail/rational_horner1_13.hpp>
  117. #include <boost/math/tools/detail/rational_horner1_14.hpp>
  118. #include <boost/math/tools/detail/rational_horner1_15.hpp>
  119. #include <boost/math/tools/detail/rational_horner1_16.hpp>
  120. #include <boost/math/tools/detail/rational_horner1_17.hpp>
  121. #include <boost/math/tools/detail/rational_horner1_18.hpp>
  122. #include <boost/math/tools/detail/rational_horner1_19.hpp>
  123. #include <boost/math/tools/detail/rational_horner1_20.hpp>
  124. #include <boost/math/tools/detail/rational_horner2_2.hpp>
  125. #include <boost/math/tools/detail/rational_horner2_3.hpp>
  126. #include <boost/math/tools/detail/rational_horner2_4.hpp>
  127. #include <boost/math/tools/detail/rational_horner2_5.hpp>
  128. #include <boost/math/tools/detail/rational_horner2_6.hpp>
  129. #include <boost/math/tools/detail/rational_horner2_7.hpp>
  130. #include <boost/math/tools/detail/rational_horner2_8.hpp>
  131. #include <boost/math/tools/detail/rational_horner2_9.hpp>
  132. #include <boost/math/tools/detail/rational_horner2_10.hpp>
  133. #include <boost/math/tools/detail/rational_horner2_11.hpp>
  134. #include <boost/math/tools/detail/rational_horner2_12.hpp>
  135. #include <boost/math/tools/detail/rational_horner2_13.hpp>
  136. #include <boost/math/tools/detail/rational_horner2_14.hpp>
  137. #include <boost/math/tools/detail/rational_horner2_15.hpp>
  138. #include <boost/math/tools/detail/rational_horner2_16.hpp>
  139. #include <boost/math/tools/detail/rational_horner2_17.hpp>
  140. #include <boost/math/tools/detail/rational_horner2_18.hpp>
  141. #include <boost/math/tools/detail/rational_horner2_19.hpp>
  142. #include <boost/math/tools/detail/rational_horner2_20.hpp>
  143. #include <boost/math/tools/detail/rational_horner3_2.hpp>
  144. #include <boost/math/tools/detail/rational_horner3_3.hpp>
  145. #include <boost/math/tools/detail/rational_horner3_4.hpp>
  146. #include <boost/math/tools/detail/rational_horner3_5.hpp>
  147. #include <boost/math/tools/detail/rational_horner3_6.hpp>
  148. #include <boost/math/tools/detail/rational_horner3_7.hpp>
  149. #include <boost/math/tools/detail/rational_horner3_8.hpp>
  150. #include <boost/math/tools/detail/rational_horner3_9.hpp>
  151. #include <boost/math/tools/detail/rational_horner3_10.hpp>
  152. #include <boost/math/tools/detail/rational_horner3_11.hpp>
  153. #include <boost/math/tools/detail/rational_horner3_12.hpp>
  154. #include <boost/math/tools/detail/rational_horner3_13.hpp>
  155. #include <boost/math/tools/detail/rational_horner3_14.hpp>
  156. #include <boost/math/tools/detail/rational_horner3_15.hpp>
  157. #include <boost/math/tools/detail/rational_horner3_16.hpp>
  158. #include <boost/math/tools/detail/rational_horner3_17.hpp>
  159. #include <boost/math/tools/detail/rational_horner3_18.hpp>
  160. #include <boost/math/tools/detail/rational_horner3_19.hpp>
  161. #include <boost/math/tools/detail/rational_horner3_20.hpp>
  162. #endif
  163. namespace boost{ namespace math{ namespace tools{
  164. //
  165. // Forward declaration to keep two phase lookup happy:
  166. //
  167. template <class T, class U>
  168. BOOST_MATH_GPU_ENABLED U evaluate_polynomial(const T* poly, U const& z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U);
  169. namespace detail{
  170. template <class T, class V, class Tag>
  171. 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)
  172. {
  173. return evaluate_polynomial(a, val, Tag::value);
  174. }
  175. } // namespace detail
  176. //
  177. // Polynomial evaluation with runtime size.
  178. // This requires a for-loop which may be more expensive than
  179. // the loop expanded versions above:
  180. //
  181. template <class T, class U>
  182. BOOST_MATH_GPU_ENABLED inline U evaluate_polynomial(const T* poly, U const& z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
  183. {
  184. BOOST_MATH_ASSERT(count > 0);
  185. U sum = static_cast<U>(poly[count - 1]);
  186. for(int i = static_cast<int>(count) - 2; i >= 0; --i)
  187. {
  188. sum *= z;
  189. sum += static_cast<U>(poly[i]);
  190. }
  191. return sum;
  192. }
  193. //
  194. // Compile time sized polynomials, just inline forwarders to the
  195. // implementations above:
  196. //
  197. template <boost::math::size_t N, class T, class V>
  198. BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_polynomial(const T(&a)[N], const V& val) BOOST_MATH_NOEXCEPT(V)
  199. {
  200. typedef boost::math::integral_constant<int, static_cast<int>(N)> tag_type;
  201. return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(nullptr));
  202. }
  203. #ifndef BOOST_MATH_HAS_NVRTC
  204. template <boost::math::size_t N, class T, class V>
  205. 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)
  206. {
  207. typedef boost::math::integral_constant<int, static_cast<int>(N)> tag_type;
  208. return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()), val, static_cast<tag_type const*>(nullptr));
  209. }
  210. #endif
  211. //
  212. // Even polynomials are trivial: just square the argument!
  213. //
  214. template <class T, class U>
  215. BOOST_MATH_GPU_ENABLED inline U evaluate_even_polynomial(const T* poly, U z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
  216. {
  217. return evaluate_polynomial(poly, U(z*z), count);
  218. }
  219. template <boost::math::size_t N, class T, class V>
  220. BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_even_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V)
  221. {
  222. return evaluate_polynomial(a, V(z*z));
  223. }
  224. #ifndef BOOST_MATH_HAS_NVRTC
  225. template <boost::math::size_t N, class T, class V>
  226. 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)
  227. {
  228. return evaluate_polynomial(a, V(z*z));
  229. }
  230. #endif
  231. //
  232. // Odd polynomials come next:
  233. //
  234. template <class T, class U>
  235. BOOST_MATH_GPU_ENABLED inline U evaluate_odd_polynomial(const T* poly, U z, boost::math::size_t count) BOOST_MATH_NOEXCEPT(U)
  236. {
  237. return poly[0] + z * evaluate_polynomial(poly+1, U(z*z), count-1);
  238. }
  239. template <boost::math::size_t N, class T, class V>
  240. BOOST_MATH_GPU_ENABLED BOOST_MATH_GPU_ENABLED inline V evaluate_odd_polynomial(const T(&a)[N], const V& z) BOOST_MATH_NOEXCEPT(V)
  241. {
  242. typedef boost::math::integral_constant<int, static_cast<int>(N-1)> tag_type;
  243. return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, V(z*z), static_cast<tag_type const*>(nullptr));
  244. }
  245. #ifndef BOOST_MATH_HAS_NVRTC
  246. template <boost::math::size_t N, class T, class V>
  247. 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)
  248. {
  249. typedef boost::math::integral_constant<int, static_cast<int>(N-1)> tag_type;
  250. 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));
  251. }
  252. #endif
  253. template <class T, class U, class V>
  254. 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);
  255. namespace detail{
  256. template <class T, class U, class V, class Tag>
  257. 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)
  258. {
  259. return boost::math::tools::evaluate_rational(num, denom, z, Tag::value);
  260. }
  261. }
  262. //
  263. // Rational functions: numerator and denominator must be
  264. // equal in size. These always have a for-loop and so may be less
  265. // efficient than evaluating a pair of polynomials. However, there
  266. // are some tricks we can use to prevent overflow that might otherwise
  267. // occur in polynomial evaluation, if z is large. This is important
  268. // in our Lanczos code for example.
  269. //
  270. template <class T, class U, class V>
  271. 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)
  272. {
  273. V z(z_);
  274. V s1, s2;
  275. if(z <= 1)
  276. {
  277. s1 = static_cast<V>(num[count-1]);
  278. s2 = static_cast<V>(denom[count-1]);
  279. for(int i = (int)count - 2; i >= 0; --i)
  280. {
  281. s1 *= z;
  282. s2 *= z;
  283. s1 += num[i];
  284. s2 += denom[i];
  285. }
  286. }
  287. else
  288. {
  289. z = 1 / z;
  290. s1 = static_cast<V>(num[0]);
  291. s2 = static_cast<V>(denom[0]);
  292. for(unsigned i = 1; i < count; ++i)
  293. {
  294. s1 *= z;
  295. s2 *= z;
  296. s1 += num[i];
  297. s2 += denom[i];
  298. }
  299. }
  300. return s1 / s2;
  301. }
  302. template <boost::math::size_t N, class T, class U, class V>
  303. 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)
  304. {
  305. return detail::evaluate_rational_c_imp(a, b, z, static_cast<const boost::math::integral_constant<int, static_cast<int>(N)>*>(nullptr));
  306. }
  307. #ifndef BOOST_MATH_HAS_NVRTC
  308. template <boost::math::size_t N, class T, class U, class V>
  309. 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)
  310. {
  311. return detail::evaluate_rational_c_imp(a.data(), b.data(), z, static_cast<boost::math::integral_constant<int, static_cast<int>(N)>*>(nullptr));
  312. }
  313. #endif
  314. } // namespace tools
  315. } // namespace math
  316. } // namespace boost
  317. #endif // BOOST_MATH_TOOLS_RATIONAL_HPP