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