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
13 #include <boost/math/tools/precision.hpp>
14 #include <boost/math/tools/complex.hpp>
15 #include <type_traits>
16 #include <cstdint>
17 #include <cmath>
18
19 namespace boost{ namespace math{ namespace tools{
20
21 namespace detail
22 {
23
24 template <typename T>
25 struct is_pair : public std::false_type{};
26
27 template <typename T, typename U>
28 struct is_pair<std::pair<T,U>> : public std::true_type{};
29
30 template <typename Gen>
31 struct fraction_traits_simple
32 {
33 using result_type = typename Gen::result_type;
34 using value_type = typename Gen::result_type;
35
aboost::math::tools::detail::fraction_traits_simple36 static result_type a(const value_type&) BOOST_MATH_NOEXCEPT(value_type)
37 {
38 return 1;
39 }
bboost::math::tools::detail::fraction_traits_simple40 static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
41 {
42 return v;
43 }
44 };
45
46 template <typename Gen>
47 struct fraction_traits_pair
48 {
49 using value_type = typename Gen::result_type;
50 using result_type = typename value_type::first_type;
51
aboost::math::tools::detail::fraction_traits_pair52 static result_type a(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
53 {
54 return v.first;
55 }
bboost::math::tools::detail::fraction_traits_pair56 static result_type b(const value_type& v) BOOST_MATH_NOEXCEPT(value_type)
57 {
58 return v.second;
59 }
60 };
61
62 template <typename Gen>
63 struct fraction_traits
64 : public std::conditional<
65 is_pair<typename Gen::result_type>::value,
66 fraction_traits_pair<Gen>,
67 fraction_traits_simple<Gen>>::type
68 {
69 };
70
71 template <typename T, bool = is_complex_type<T>::value>
72 struct tiny_value
73 {
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.
getboost::math::tools::detail::tiny_value77 static T get() {
78 return 16*tools::min_value<T>();
79 }
80 };
81 template <typename T>
82 struct tiny_value<T, true>
83 {
84 using value_type = typename T::value_type;
getboost::math::tools::detail::tiny_value85 static T get() {
86 return 16*tools::min_value<value_type>();
87 }
88 };
89
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 //
104 // Note that the first a0 returned by generator Gen is discarded.
105 //
106 template <typename Gen, typename U>
continued_fraction_b(Gen & g,const U & factor,std::uintmax_t & max_terms)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>()()))
109 {
110 BOOST_MATH_STD_USING // ADL of std names
111
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;
117
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);
122
123 value_type v = g();
124
125 result_type f, C, D, delta;
126 f = traits::b(v);
127 if(f == zero)
128 f = tiny;
129 C = f;
130 D = 0;
131
132 std::uintmax_t counter(max_terms);
133 do{
134 v = g();
135 D = traits::b(v) + traits::a(v) * D;
136 if(D == result_type(0))
137 D = tiny;
138 C = traits::b(v) + traits::a(v) / C;
139 if(C == zero)
140 C = tiny;
141 D = one/D;
142 delta = C*D;
143 f = f * delta;
144 }while((abs(delta - one) > terminator) && --counter);
145
146 max_terms = max_terms - counter;
147
148 return f;
149 }
150
151 template <typename Gen, typename U>
continued_fraction_b(Gen & g,const U & factor)152 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
153 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
154 {
155 std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
156 return continued_fraction_b(g, factor, max_terms);
157 }
158
159 template <typename Gen>
continued_fraction_b(Gen & g,int bits)160 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
161 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
162 {
163 BOOST_MATH_STD_USING // ADL of std names
164
165 using traits = detail::fraction_traits<Gen>;
166 using result_type = typename traits::result_type;
167
168 result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
169 std::uintmax_t max_terms = (std::numeric_limits<std::uintmax_t>::max)();
170 return continued_fraction_b(g, factor, max_terms);
171 }
172
173 template <typename Gen>
continued_fraction_b(Gen & g,int bits,std::uintmax_t & max_terms)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>()()))
176 {
177 BOOST_MATH_STD_USING // ADL of std names
178
179 using traits = detail::fraction_traits<Gen>;
180 using result_type = typename traits::result_type;
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 //
200 template <typename Gen, typename U>
continued_fraction_a(Gen & g,const U & factor,std::uintmax_t & max_terms)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>()()))
203 {
204 BOOST_MATH_STD_USING // ADL of std names
205
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;
211
212 integer_type const zero(0), one(1);
213
214 result_type tiny = detail::tiny_value<result_type>::get();
215 scalar_type terminator = abs(factor);
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);
222 if(f == zero)
223 f = tiny;
224 C = f;
225 D = 0;
226
227 std::uintmax_t counter(max_terms);
228
229 do{
230 v = g();
231 D = traits::b(v) + traits::a(v) * D;
232 if(D == zero)
233 D = tiny;
234 C = traits::b(v) + traits::a(v) / C;
235 if(C == zero)
236 C = tiny;
237 D = one/D;
238 delta = C*D;
239 f = f * delta;
240 }while((abs(delta - one) > terminator) && --counter);
241
242 max_terms = max_terms - counter;
243
244 return a0/f;
245 }
246
247 template <typename Gen, typename U>
continued_fraction_a(Gen & g,const U & factor)248 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
249 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
250 {
251 std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
252 return continued_fraction_a(g, factor, max_iter);
253 }
254
255 template <typename Gen>
continued_fraction_a(Gen & g,int bits)256 inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
257 noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()()))
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);
265 std::uintmax_t max_iter = (std::numeric_limits<std::uintmax_t>::max)();
266
267 return continued_fraction_a(g, factor, max_iter);
268 }
269
270 template <typename Gen>
continued_fraction_a(Gen & g,int bits,std::uintmax_t & max_terms)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>()()))
273 {
274 BOOST_MATH_STD_USING // ADL of std names
275
276 using traits = detail::fraction_traits<Gen>;
277 using result_type = typename traits::result_type;
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
288