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