fraction.hpp 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252
  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&)
  28. {
  29. return 1;
  30. }
  31. static result_type b(const value_type& v)
  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)
  42. {
  43. return v.first;
  44. }
  45. static result_type b(const value_type& v)
  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. {
  76. BOOST_MATH_STD_USING // ADL of std names
  77. typedef detail::fraction_traits<Gen> traits;
  78. typedef typename traits::result_type result_type;
  79. typedef typename traits::value_type value_type;
  80. result_type tiny = tools::min_value<result_type>();
  81. value_type v = g();
  82. result_type f, C, D, delta;
  83. f = traits::b(v);
  84. if(f == 0)
  85. f = tiny;
  86. C = f;
  87. D = 0;
  88. boost::uintmax_t counter(max_terms);
  89. do{
  90. v = g();
  91. D = traits::b(v) + traits::a(v) * D;
  92. if(D == 0)
  93. D = tiny;
  94. C = traits::b(v) + traits::a(v) / C;
  95. if(C == 0)
  96. C = tiny;
  97. D = 1/D;
  98. delta = C*D;
  99. f = f * delta;
  100. }while((fabs(delta - 1) > factor) && --counter);
  101. max_terms = max_terms - counter;
  102. return f;
  103. }
  104. template <class Gen, class U>
  105. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
  106. {
  107. boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
  108. return continued_fraction_b(g, factor, max_terms);
  109. }
  110. template <class Gen>
  111. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
  112. {
  113. BOOST_MATH_STD_USING // ADL of std names
  114. typedef detail::fraction_traits<Gen> traits;
  115. typedef typename traits::result_type result_type;
  116. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  117. boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)();
  118. return continued_fraction_b(g, factor, max_terms);
  119. }
  120. template <class Gen>
  121. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
  122. {
  123. BOOST_MATH_STD_USING // ADL of std names
  124. typedef detail::fraction_traits<Gen> traits;
  125. typedef typename traits::result_type result_type;
  126. result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
  127. return continued_fraction_b(g, factor, max_terms);
  128. }
  129. //
  130. // continued_fraction_a
  131. // Evaluates:
  132. //
  133. // a1
  134. // ---------------
  135. // b1 + a2
  136. // ----------
  137. // b2 + a3
  138. // -----
  139. // b3 + ...
  140. //
  141. // Note that the first a1 and b1 returned by generator Gen are both used.
  142. //
  143. template <class Gen, class U>
  144. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms)
  145. {
  146. BOOST_MATH_STD_USING // ADL of std names
  147. typedef detail::fraction_traits<Gen> traits;
  148. typedef typename traits::result_type result_type;
  149. typedef typename traits::value_type value_type;
  150. result_type tiny = tools::min_value<result_type>();
  151. value_type v = g();
  152. result_type f, C, D, delta, a0;
  153. f = traits::b(v);
  154. a0 = traits::a(v);
  155. if(f == 0)
  156. f = tiny;
  157. C = f;
  158. D = 0;
  159. boost::uintmax_t counter(max_terms);
  160. do{
  161. v = g();
  162. D = traits::b(v) + traits::a(v) * D;
  163. if(D == 0)
  164. D = tiny;
  165. C = traits::b(v) + traits::a(v) / C;
  166. if(C == 0)
  167. C = tiny;
  168. D = 1/D;
  169. delta = C*D;
  170. f = f * delta;
  171. }while((fabs(delta - 1) > factor) && --counter);
  172. max_terms = max_terms - counter;
  173. return a0/f;
  174. }
  175. template <class Gen, class U>
  176. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
  177. {
  178. boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
  179. return continued_fraction_a(g, factor, max_iter);
  180. }
  181. template <class Gen>
  182. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
  183. {
  184. BOOST_MATH_STD_USING // ADL of std names
  185. typedef detail::fraction_traits<Gen> traits;
  186. typedef typename traits::result_type result_type;
  187. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  188. boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)();
  189. return continued_fraction_a(g, factor, max_iter);
  190. }
  191. template <class Gen>
  192. inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
  193. {
  194. BOOST_MATH_STD_USING // ADL of std names
  195. typedef detail::fraction_traits<Gen> traits;
  196. typedef typename traits::result_type result_type;
  197. result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
  198. return continued_fraction_a(g, factor, max_terms);
  199. }
  200. } // namespace tools
  201. } // namespace math
  202. } // namespace boost
  203. #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED