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