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