fraction.hpp 8.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287
  1. // (C) Copyright John Maddock 2005-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_FRACTION_INCLUDED
  6. #define BOOST_MATH_TOOLS_FRACTION_INCLUDED
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/math/tools/precision.hpp>
  11. #include <boost/math/tools/complex.hpp>
  12. #include <type_traits>
  13. #include <cstdint>
  14. #include <cmath>
  15. namespace boost{ namespace math{ namespace tools{
  16. namespace detail
  17. {
  18. template <typename T>
  19. struct is_pair : public std::false_type{};
  20. template <typename T, typename U>
  21. struct is_pair<std::pair<T,U>> : public std::true_type{};
  22. template <typename Gen>
  23. struct fraction_traits_simple
  24. {
  25. using result_type = typename Gen::result_type;
  26. using value_type = typename Gen::result_type;
  27. static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
  28. {
  29. return 1;
  30. }
  31. static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  32. {
  33. return v;
  34. }
  35. };
  36. template <typename Gen>
  37. struct fraction_traits_pair
  38. {
  39. using value_type = typename Gen::result_type;
  40. using result_type = typename value_type::first_type;
  41. static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  42. {
  43. return v.first;
  44. }
  45. static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
  46. {
  47. return v.second;
  48. }
  49. };
  50. template <typename Gen>
  51. struct fraction_traits
  52. : public std::conditional<
  53. is_pair<typename Gen::result_type>::value,
  54. fraction_traits_pair<Gen>,
  55. fraction_traits_simple<Gen>>::type
  56. {
  57. };
  58. template <typename T, bool = is_complex_type<T>::value>
  59. struct tiny_value
  60. {
  61. // For float, double, and long double, 1/min_value<T>() is finite.
  62. // But for mpfr_float and cpp_bin_float, 1/min_value<T>() is inf.
  63. // Multiply the min by 16 so that the reciprocal doesn't overflow.
  64. static T get() {
  65. return 16*tools::min_value<T>();
  66. }
  67. };
  68. template <typename T>
  69. struct tiny_value<T, true>
  70. {
  71. using value_type = typename T::value_type;
  72. static T get() {
  73. return 16*tools::min_value<value_type>();
  74. }
  75. };
  76. } // namespace detail
  77. //
  78. // continued_fraction_b
  79. // Evaluates:
  80. //
  81. // b0 + a1
  82. // ---------------
  83. // b1 + a2
  84. // ----------
  85. // b2 + a3
  86. // -----
  87. // b3 + ...
  88. //
  89. // Note that the first a0 returned by generator Gen is discarded.
  90. //
  91. template <typename Gen, typename U>
  92. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, std::uintmax_t& max_terms)
  93. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  94. {
  95. BOOST_MATH_STD_USING // ADL of std names
  96. using traits = detail::fraction_traits<Gen>;
  97. using result_type = typename traits::result_type;
  98. using value_type = typename traits::value_type;
  99. using integer_type = typename integer_scalar_type<result_type>::type;
  100. using scalar_type = typename scalar_type<result_type>::type;
  101. integer_type const zero(0), one(1);
  102. result_type tiny = detail::tiny_value<result_type>::get();
  103. scalar_type terminator = abs(factor);
  104. value_type v = g();
  105. result_type f, C, D, delta;
  106. f = traits::b(v);
  107. if(f == zero)
  108. f = tiny;
  109. C = f;
  110. D = 0;
  111. std::uintmax_t counter(max_terms);
  112. do{
  113. v = g();
  114. D = traits::b(v) + traits::a(v) * D;
  115. if(D == result_type(0))
  116. D = tiny;
  117. C = traits::b(v) + traits::a(v) / C;
  118. if(C == zero)
  119. C = tiny;
  120. D = one/D;
  121. delta = C*D;
  122. f = f * delta;
  123. }while((abs(delta - one) > terminator) && --counter);
  124. max_terms = max_terms - counter;
  125. return f;
  126. }
  127. template <typename Gen, typename U>
  128. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
  129. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  130. {
  131. std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
  132. return continued_fraction_b(g, factor, max_terms);
  133. }
  134. template <typename Gen>
  135. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
  136. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  137. {
  138. BOOST_MATH_STD_USING // ADL of std names
  139. using traits = detail::fraction_traits<Gen>;
  140. using result_type = typename traits::result_type;
  141. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  142. std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
  143. return continued_fraction_b(g, factor, max_terms);
  144. }
  145. template <typename Gen>
  146. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, std::uintmax_t& max_terms)
  147. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  148. {
  149. BOOST_MATH_STD_USING // ADL of std names
  150. using traits = detail::fraction_traits<Gen>;
  151. using result_type = typename traits::result_type;
  152. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  153. return continued_fraction_b(g, factor, max_terms);
  154. }
  155. //
  156. // continued_fraction_a
  157. // Evaluates:
  158. //
  159. // a1
  160. // ---------------
  161. // b1 + a2
  162. // ----------
  163. // b2 + a3
  164. // -----
  165. // b3 + ...
  166. //
  167. // Note that the first a1 and b1 returned by generator Gen are both used.
  168. //
  169. template <typename Gen, typename U>
  170. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, std::uintmax_t& max_terms)
  171. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  172. {
  173. BOOST_MATH_STD_USING // ADL of std names
  174. using traits = detail::fraction_traits<Gen>;
  175. using result_type = typename traits::result_type;
  176. using value_type = typename traits::value_type;
  177. using integer_type = typename integer_scalar_type<result_type>::type;
  178. using scalar_type = typename scalar_type<result_type>::type;
  179. integer_type const zero(0), one(1);
  180. result_type tiny = detail::tiny_value<result_type>::get();
  181. scalar_type terminator = abs(factor);
  182. value_type v = g();
  183. result_type f, C, D, delta, a0;
  184. f = traits::b(v);
  185. a0 = traits::a(v);
  186. if(f == zero)
  187. f = tiny;
  188. C = f;
  189. D = 0;
  190. std::uintmax_t counter(max_terms);
  191. do{
  192. v = g();
  193. D = traits::b(v) + traits::a(v) * D;
  194. if(D == zero)
  195. D = tiny;
  196. C = traits::b(v) + traits::a(v) / C;
  197. if(C == zero)
  198. C = tiny;
  199. D = one/D;
  200. delta = C*D;
  201. f = f * delta;
  202. }while((abs(delta - one) > terminator) && --counter);
  203. max_terms = max_terms - counter;
  204. return a0/f;
  205. }
  206. template <typename Gen, typename U>
  207. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
  208. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  209. {
  210. std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
  211. return continued_fraction_a(g, factor, max_iter);
  212. }
  213. template <typename Gen>
  214. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
  215. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  216. {
  217. BOOST_MATH_STD_USING // ADL of std names
  218. typedef detail::fraction_traits<Gen> traits;
  219. typedef typename traits::result_type result_type;
  220. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  221. std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
  222. return continued_fraction_a(g, factor, max_iter);
  223. }
  224. template <typename Gen>
  225. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, std::uintmax_t& max_terms)
  226. noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  227. {
  228. BOOST_MATH_STD_USING // ADL of std names
  229. using traits = detail::fraction_traits<Gen>;
  230. using result_type = typename traits::result_type;
  231. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  232. return continued_fraction_a(g, factor, max_terms);
  233. }
  234. } // namespace tools
  235. } // namespace math
  236. } // namespace boost
  237. #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED