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