fraction.hpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355
  1. // (C) Copyright John Maddock 2005-2006.
  2. // (C) Copyright Matt Borland 2024.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_MATH_TOOLS_FRACTION_INCLUDED
  7. #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
  8. #ifdef _MSC_VER
  9. #pragma once
  10. #endif
  11. #include <boost/math/tools/config.hpp>
  12. #include <boost/math/tools/type_traits.hpp>
  13. #include <boost/math/tools/numeric_limits.hpp>
  14. #include <boost/math/tools/tuple.hpp>
  15. #include <boost/math/tools/precision.hpp>
  16. #include <boost/math/tools/complex.hpp>
  17. #include <boost/math/tools/cstdint.hpp>
  18. namespace boost{ namespace math{ namespace tools{
  19. namespace detail
  20. {
  21. template <typename T>
  22. struct is_pair : public boost::math::false_type{};
  23. template <typename T, typename U>
  24. struct is_pair<boost::math::pair<T,U>> : public boost::math::true_type{};
  25. template <typename Gen>
  26. struct fraction_traits_simple
  27. {
  28. using result_type = typename Gen::result_type;
  29. using value_type = typename Gen::result_type;
  30. BOOST_MATH_GPU_ENABLED static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
  31. {
  32. return 1;
  33. }
  34. BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  35. {
  36. return v;
  37. }
  38. };
  39. template <typename Gen>
  40. struct fraction_traits_pair
  41. {
  42. using value_type = typename Gen::result_type;
  43. using result_type = typename value_type::first_type;
  44. BOOST_MATH_GPU_ENABLED static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  45. {
  46. return v.first;
  47. }
  48. BOOST_MATH_GPU_ENABLED static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  49. {
  50. return v.second;
  51. }
  52. };
  53. template <typename Gen>
  54. struct fraction_traits
  55. : public boost::math::conditional<
  56. is_pair<typename Gen::result_type>::value,
  57. fraction_traits_pair<Gen>,
  58. fraction_traits_simple<Gen>>::type
  59. {
  60. };
  61. template <typename T, bool = is_complex_type<T>::value>
  62. struct tiny_value
  63. {
  64. // For float, double, and long double, 1/min_value<T>() is finite.
  65. // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.
  66. // Multiply the min by 16 so that the reciprocal doesn't overflow.
  67. BOOST_MATH_GPU_ENABLED static T get() {
  68. return 16*tools::min_value<T>();
  69. }
  70. };
  71. template <typename T>
  72. struct tiny_value<T, true>
  73. {
  74. using value_type = typename T::value_type;
  75. BOOST_MATH_GPU_ENABLED static T get() {
  76. return 16*tools::min_value<value_type>();
  77. }
  78. };
  79. } // namespace detail
  80. namespace detail {
  81. //
  82. // continued_fraction_b
  83. // Evaluates:
  84. //
  85. // b0 + a1
  86. // ---------------
  87. // b1 + a2
  88. // ----------
  89. // b2 + a3
  90. // -----
  91. // b3 + ...
  92. //
  93. // Note that the first a0 returned by generator Gen is discarded.
  94. //
  95. template <typename Gen, typename U>
  96. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
  97. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  98. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  99. // SYCL can not handle this condition so we only check float on that platform
  100. && noexcept(std::declval<Gen>()())
  101. #endif
  102. )
  103. {
  104. BOOST_MATH_STD_USING // ADL of std names
  105. using traits = detail::fraction_traits<Gen>;
  106. using result_type = typename traits::result_type;
  107. using value_type = typename traits::value_type;
  108. using integer_type = typename integer_scalar_type<result_type>::type;
  109. using scalar_type = typename scalar_type<result_type>::type;
  110. integer_type const zero(0), one(1);
  111. result_type tiny = detail::tiny_value<result_type>::get();
  112. scalar_type terminator = abs(factor);
  113. value_type v = g();
  114. result_type f, C, D, delta;
  115. f = traits::b(v);
  116. if(f == zero)
  117. f = tiny;
  118. C = f;
  119. D = 0;
  120. boost::math::uintmax_t counter(max_terms);
  121. do{
  122. v = g();
  123. D = traits::b(v) + traits::a(v) * D;
  124. if(D == result_type(0))
  125. D = tiny;
  126. C = traits::b(v) + traits::a(v) / C;
  127. if(C == zero)
  128. C = tiny;
  129. D = one/D;
  130. delta = C*D;
  131. f = f * delta;
  132. }while((abs(delta - one) > terminator) && --counter);
  133. max_terms = max_terms - counter;
  134. return f;
  135. }
  136. } // namespace detail
  137. template <typename Gen, typename U>
  138. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
  139. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  140. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  141. && noexcept(std::declval<Gen>()())
  142. #endif
  143. )
  144. {
  145. return detail::continued_fraction_b_impl(g, factor, max_terms);
  146. }
  147. template <typename Gen, typename U>
  148. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
  149. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  150. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  151. && noexcept(std::declval<Gen>()())
  152. #endif
  153. )
  154. {
  155. boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  156. return detail::continued_fraction_b_impl(g, factor, max_terms);
  157. }
  158. template <typename Gen>
  159. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
  160. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  161. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  162. && noexcept(std::declval<Gen>()())
  163. #endif
  164. )
  165. {
  166. BOOST_MATH_STD_USING // ADL of std names
  167. using traits = detail::fraction_traits<Gen>;
  168. using result_type = typename traits::result_type;
  169. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  170. boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  171. return detail::continued_fraction_b_impl(g, factor, max_terms);
  172. }
  173. template <typename Gen>
  174. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms)
  175. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  176. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  177. && noexcept(std::declval<Gen>()())
  178. #endif
  179. )
  180. {
  181. BOOST_MATH_STD_USING // ADL of std names
  182. using traits = detail::fraction_traits<Gen>;
  183. using result_type = typename traits::result_type;
  184. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  185. return detail::continued_fraction_b_impl(g, factor, max_terms);
  186. }
  187. namespace detail {
  188. //
  189. // continued_fraction_a
  190. // Evaluates:
  191. //
  192. // a1
  193. // ---------------
  194. // b1 + a2
  195. // ----------
  196. // b2 + a3
  197. // -----
  198. // b3 + ...
  199. //
  200. // Note that the first a1 and b1 returned by generator Gen are both used.
  201. //
  202. template <typename Gen, typename U>
  203. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
  204. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  205. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  206. && noexcept(std::declval<Gen>()())
  207. #endif
  208. )
  209. {
  210. BOOST_MATH_STD_USING // ADL of std names
  211. using traits = detail::fraction_traits<Gen>;
  212. using result_type = typename traits::result_type;
  213. using value_type = typename traits::value_type;
  214. using integer_type = typename integer_scalar_type<result_type>::type;
  215. using scalar_type = typename scalar_type<result_type>::type;
  216. integer_type const zero(0), one(1);
  217. result_type tiny = detail::tiny_value<result_type>::get();
  218. scalar_type terminator = abs(factor);
  219. value_type v = g();
  220. result_type f, C, D, delta, a0;
  221. f = traits::b(v);
  222. a0 = traits::a(v);
  223. if(f == zero)
  224. f = tiny;
  225. C = f;
  226. D = 0;
  227. boost::math::uintmax_t counter(max_terms);
  228. do{
  229. v = g();
  230. D = traits::b(v) + traits::a(v) * D;
  231. if(D == zero)
  232. D = tiny;
  233. C = traits::b(v) + traits::a(v) / C;
  234. if(C == zero)
  235. C = tiny;
  236. D = one/D;
  237. delta = C*D;
  238. f = f * delta;
  239. }while((abs(delta - one) > terminator) && --counter);
  240. max_terms = max_terms - counter;
  241. return a0/f;
  242. }
  243. } // namespace detail
  244. template <typename Gen, typename U>
  245. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
  246. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  247. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  248. && noexcept(std::declval<Gen>()())
  249. #endif
  250. )
  251. {
  252. return detail::continued_fraction_a_impl(g, factor, max_terms);
  253. }
  254. template <typename Gen, typename U>
  255. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
  256. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  257. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  258. && noexcept(std::declval<Gen>()())
  259. #endif
  260. )
  261. {
  262. boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  263. return detail::continued_fraction_a_impl(g, factor, max_iter);
  264. }
  265. template <typename Gen>
  266. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
  267. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  268. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  269. && noexcept(std::declval<Gen>()())
  270. #endif
  271. )
  272. {
  273. BOOST_MATH_STD_USING // ADL of std names
  274. typedef detail::fraction_traits<Gen> traits;
  275. typedef typename traits::result_type result_type;
  276. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  277. boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  278. return detail::continued_fraction_a_impl(g, factor, max_iter);
  279. }
  280. template <typename Gen>
  281. BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms)
  282. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
  283. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  284. && noexcept(std::declval<Gen>()())
  285. #endif
  286. )
  287. {
  288. BOOST_MATH_STD_USING // ADL of std names
  289. using traits = detail::fraction_traits<Gen>;
  290. using result_type = typename traits::result_type;
  291. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  292. return detail::continued_fraction_a_impl(g, factor, max_terms);
  293. }
  294. } // namespace tools
  295. } // namespace math
  296. } // namespace boost
  297. #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED