]>
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 | ||
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 | ||
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 | } | |
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 | ||
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 | } | |
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> | |
88 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::uintmax_t& max_terms) | |
89 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
90 | { | |
91 | BOOST_MATH_STD_USING // ADL of std names | |
92 | ||
93 | typedef detail::fraction_traits<Gen> traits; | |
94 | typedef typename traits::result_type result_type; | |
95 | typedef typename traits::value_type value_type; | |
96 | ||
97 | result_type tiny = tools::min_value<result_type>(); | |
98 | ||
99 | value_type v = g(); | |
100 | ||
101 | result_type f, C, D, delta; | |
102 | f = traits::b(v); | |
103 | if(f == 0) | |
104 | f = tiny; | |
105 | C = f; | |
106 | D = 0; | |
107 | ||
108 | boost::uintmax_t counter(max_terms); | |
109 | ||
110 | do{ | |
111 | v = g(); | |
112 | D = traits::b(v) + traits::a(v) * D; | |
113 | if(D == 0) | |
114 | D = tiny; | |
115 | C = traits::b(v) + traits::a(v) / C; | |
116 | if(C == 0) | |
117 | C = tiny; | |
118 | D = 1/D; | |
119 | delta = C*D; | |
120 | f = f * delta; | |
121 | }while((fabs(delta - 1) > factor) && --counter); | |
122 | ||
123 | max_terms = max_terms - counter; | |
124 | ||
125 | return f; | |
126 | } | |
127 | ||
128 | template <class Gen, class U> | |
129 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor) | |
130 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
131 | { | |
132 | boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)(); | |
133 | return continued_fraction_b(g, factor, max_terms); | |
134 | } | |
135 | ||
136 | template <class Gen> | |
137 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits) | |
138 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
139 | { | |
140 | BOOST_MATH_STD_USING // ADL of std names | |
141 | ||
142 | typedef detail::fraction_traits<Gen> traits; | |
143 | typedef typename traits::result_type result_type; | |
144 | ||
145 | result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits); | |
146 | boost::uintmax_t max_terms = (std::numeric_limits<boost::uintmax_t>::max)(); | |
147 | return continued_fraction_b(g, factor, max_terms); | |
148 | } | |
149 | ||
150 | template <class Gen> | |
151 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms) | |
152 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
153 | { | |
154 | BOOST_MATH_STD_USING // ADL of std names | |
155 | ||
156 | typedef detail::fraction_traits<Gen> traits; | |
157 | typedef typename traits::result_type result_type; | |
158 | ||
159 | result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits); | |
160 | return continued_fraction_b(g, factor, max_terms); | |
161 | } | |
162 | ||
163 | // | |
164 | // continued_fraction_a | |
165 | // Evaluates: | |
166 | // | |
167 | // a1 | |
168 | // --------------- | |
169 | // b1 + a2 | |
170 | // ---------- | |
171 | // b2 + a3 | |
172 | // ----- | |
173 | // b3 + ... | |
174 | // | |
175 | // Note that the first a1 and b1 returned by generator Gen are both used. | |
176 | // | |
177 | template <class Gen, class U> | |
178 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::uintmax_t& max_terms) | |
179 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
180 | { | |
181 | BOOST_MATH_STD_USING // ADL of std names | |
182 | ||
183 | typedef detail::fraction_traits<Gen> traits; | |
184 | typedef typename traits::result_type result_type; | |
185 | typedef typename traits::value_type value_type; | |
186 | ||
187 | result_type tiny = tools::min_value<result_type>(); | |
188 | ||
189 | value_type v = g(); | |
190 | ||
191 | result_type f, C, D, delta, a0; | |
192 | f = traits::b(v); | |
193 | a0 = traits::a(v); | |
194 | if(f == 0) | |
195 | f = tiny; | |
196 | C = f; | |
197 | D = 0; | |
198 | ||
199 | boost::uintmax_t counter(max_terms); | |
200 | ||
201 | do{ | |
202 | v = g(); | |
203 | D = traits::b(v) + traits::a(v) * D; | |
204 | if(D == 0) | |
205 | D = tiny; | |
206 | C = traits::b(v) + traits::a(v) / C; | |
207 | if(C == 0) | |
208 | C = tiny; | |
209 | D = 1/D; | |
210 | delta = C*D; | |
211 | f = f * delta; | |
212 | }while((fabs(delta - 1) > factor) && --counter); | |
213 | ||
214 | max_terms = max_terms - counter; | |
215 | ||
216 | return a0/f; | |
217 | } | |
218 | ||
219 | template <class Gen, class U> | |
220 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor) | |
221 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
222 | { | |
223 | boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)(); | |
224 | return continued_fraction_a(g, factor, max_iter); | |
225 | } | |
226 | ||
227 | template <class Gen> | |
228 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits) | |
229 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
230 | { | |
231 | BOOST_MATH_STD_USING // ADL of std names | |
232 | ||
233 | typedef detail::fraction_traits<Gen> traits; | |
234 | typedef typename traits::result_type result_type; | |
235 | ||
236 | result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits); | |
237 | boost::uintmax_t max_iter = (std::numeric_limits<boost::uintmax_t>::max)(); | |
238 | ||
239 | return continued_fraction_a(g, factor, max_iter); | |
240 | } | |
241 | ||
242 | template <class Gen> | |
243 | inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms) | |
244 | BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) && noexcept(std::declval<Gen>()())) | |
245 | { | |
246 | BOOST_MATH_STD_USING // ADL of std names | |
247 | ||
248 | typedef detail::fraction_traits<Gen> traits; | |
249 | typedef typename traits::result_type result_type; | |
250 | ||
251 | result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits); | |
252 | return continued_fraction_a(g, factor, max_terms); | |
253 | } | |
254 | ||
255 | } // namespace tools | |
256 | } // namespace math | |
257 | } // namespace boost | |
258 | ||
259 | #endif // BOOST_MATH_TOOLS_FRACTION_INCLUDED | |
260 |