fraction.hpp 7.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260
  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/config/no_tr1/cmath.hpp>
  11. #include <boost/cstdint.hpp>
  12. #include <boost/type_traits/integral_constant.hpp>
  13. #include <boost/mpl/if.hpp>
  14. #include <boost/math/tools/precision.hpp>
  15. namespace boost{ namespace math{ namespace tools{
  16. namespace detail
  17. {
  18. template <class T>
  19. struct is_pair : public boost::false_type{};
  20. template <class T, class U>
  21. struct is_pair<std::pair<T,U> > : public boost::true_type{};
  22. template <class Gen>
  23. struct fraction_traits_simple
  24. {
  25. typedef typename Gen::result_type result_type;
  26. typedef typename Gen::result_type value_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 <class Gen>
  37. struct fraction_traits_pair
  38. {
  39. typedef typename Gen::result_type value_type;
  40. typedef typename value_type::first_type result_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 <class Gen>
  51. struct fraction_traits
  52. : public boost::mpl::if_c<
  53. is_pair<typename Gen::result_type>::value,
  54. fraction_traits_pair<Gen>,
  55. fraction_traits_simple<Gen> >::type
  56. {
  57. };
  58. } // namespace detail
  59. //
  60. // continued_fraction_b
  61. // Evaluates:
  62. //
  63. // b0 + a1
  64. // ---------------
  65. // b1 + a2
  66. // ----------
  67. // b2 + a3
  68. // -----
  69. // b3 + ...
  70. //
  71. // Note that the first a0 returned by generator Gen is disarded.
  72. //
  73. template <class Gen, class U>
  74. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::uintmax_t& max_terms)
  75. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  76. {
  77. BOOST_MATH_STD_USING // ADL of std names
  78. typedef detail::fraction_traits<Gen> traits;
  79. typedef typename traits::result_type result_type;
  80. typedef typename traits::value_type value_type;
  81. result_type tiny = tools::min_value<result_type>();
  82. value_type v = g();
  83. result_type f, C, D, delta;
  84. f = traits::b(v);
  85. if(f == 0)
  86. f = tiny;
  87. C = f;
  88. D = 0;
  89. boost::uintmax_t counter(max_terms);
  90. do{
  91. v = g();
  92. D = traits::b(v) + traits::a(v) * D;
  93. if(D == 0)
  94. D = tiny;
  95. C = traits::b(v) + traits::a(v) / C;
  96. if(C == 0)
  97. C = tiny;
  98. D = 1/D;
  99. delta = C*D;
  100. f = f * delta;
  101. }while((fabs(delta - 1) > factor) && --counter);
  102. max_terms = max_terms - counter;
  103. return f;
  104. }
  105. template <class Gen, class U>
  106. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
  107. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  108. {
  109. boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
  110. return continued_fraction_b(g, factor, max_terms);
  111. }
  112. template <class Gen>
  113. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
  114. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  115. {
  116. BOOST_MATH_STD_USING // ADL of std names
  117. typedef detail::fraction_traits<Gen> traits;
  118. typedef typename traits::result_type result_type;
  119. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  120. boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
  121. return continued_fraction_b(g, factor, max_terms);
  122. }
  123. template <class Gen>
  124. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
  125. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  126. {
  127. BOOST_MATH_STD_USING // ADL of std names
  128. typedef detail::fraction_traits<Gen> traits;
  129. typedef typename traits::result_type result_type;
  130. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  131. return continued_fraction_b(g, factor, max_terms);
  132. }
  133. //
  134. // continued_fraction_a
  135. // Evaluates:
  136. //
  137. // a1
  138. // ---------------
  139. // b1 + a2
  140. // ----------
  141. // b2 + a3
  142. // -----
  143. // b3 + ...
  144. //
  145. // Note that the first a1 and b1 returned by generator Gen are both used.
  146. //
  147. template <class Gen, class U>
  148. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms)
  149. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  150. {
  151. BOOST_MATH_STD_USING // ADL of std names
  152. typedef detail::fraction_traits<Gen> traits;
  153. typedef typename traits::result_type result_type;
  154. typedef typename traits::value_type value_type;
  155. result_type tiny = tools::min_value<result_type>();
  156. value_type v = g();
  157. result_type f, C, D, delta, a0;
  158. f = traits::b(v);
  159. a0 = traits::a(v);
  160. if(f == 0)
  161. f = tiny;
  162. C = f;
  163. D = 0;
  164. boost::uintmax_t counter(max_terms);
  165. do{
  166. v = g();
  167. D = traits::b(v) + traits::a(v) * D;
  168. if(D == 0)
  169. D = tiny;
  170. C = traits::b(v) + traits::a(v) / C;
  171. if(C == 0)
  172. C = tiny;
  173. D = 1/D;
  174. delta = C*D;
  175. f = f * delta;
  176. }while((fabs(delta - 1) > factor) && --counter);
  177. max_terms = max_terms - counter;
  178. return a0/f;
  179. }
  180. template <class Gen, class U>
  181. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
  182. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  183. {
  184. boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
  185. return continued_fraction_a(g, factor, max_iter);
  186. }
  187. template <class Gen>
  188. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
  189. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  190. {
  191. BOOST_MATH_STD_USING // ADL of std names
  192. typedef detail::fraction_traits<Gen> traits;
  193. typedef typename traits::result_type result_type;
  194. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  195. boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
  196. return continued_fraction_a(g, factor, max_iter);
  197. }
  198. template <class Gen>
  199. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
  200. BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
  201. {
  202. BOOST_MATH_STD_USING // ADL of std names
  203. typedef detail::fraction_traits<Gen> traits;
  204. typedef typename traits::result_type result_type;
  205. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  206. return continued_fraction_a(g, factor, max_terms);
  207. }
  208. } // namespace tools
  209. } // namespace math
  210. } // namespace boost
  211. #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED