]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright Jeremy William Murphy 2016. |
2 | ||
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. (See accompanying file | |
5 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #ifndef BOOST_MATH_COMMON_FACTOR_RT_HPP | |
8 | #define BOOST_MATH_COMMON_FACTOR_RT_HPP | |
9 | ||
10 | #include <boost/assert.hpp> | |
11 | #include <boost/core/enable_if.hpp> | |
12 | #include <boost/mpl/and.hpp> | |
13 | #include <boost/type_traits.hpp> | |
14 | ||
15 | #include <boost/config.hpp> // for BOOST_NESTED_TEMPLATE, etc. | |
16 | #include <boost/limits.hpp> // for std::numeric_limits | |
17 | #include <climits> // for CHAR_MIN | |
18 | #include <boost/detail/workaround.hpp> | |
19 | #include <iterator> | |
20 | #include <algorithm> | |
21 | #include <limits> | |
22 | ||
23 | #if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64)) | |
24 | #include <intrin.h> | |
25 | #endif | |
26 | ||
27 | #ifdef BOOST_MSVC | |
28 | #pragma warning(push) | |
29 | #pragma warning(disable:4127 4244) // Conditional expression is constant | |
30 | #endif | |
31 | ||
32 | namespace boost { | |
33 | namespace math { | |
34 | ||
35 | template <class T, bool a = is_unsigned<T>::value || (std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_signed)> | |
36 | struct gcd_traits_abs_defaults | |
37 | { | |
38 | inline static const T& abs(const T& val) { return val; } | |
39 | }; | |
40 | template <class T> | |
41 | struct gcd_traits_abs_defaults<T, false> | |
42 | { | |
43 | inline static T abs(const T& val) | |
44 | { | |
45 | using std::abs; | |
46 | return abs(val); | |
47 | } | |
48 | }; | |
49 | ||
50 | template <class T> | |
51 | struct gcd_traits_defaults : public gcd_traits_abs_defaults<T> | |
52 | { | |
53 | BOOST_FORCEINLINE static unsigned make_odd(T& val) | |
54 | { | |
55 | unsigned r = 0; | |
56 | while(!(val & 1u)) | |
57 | { | |
58 | val >>= 1; | |
59 | ++r; | |
60 | } | |
61 | return r; | |
62 | } | |
63 | inline static bool less(const T& a, const T& b) | |
64 | { | |
65 | return a < b; | |
66 | } | |
67 | ||
68 | enum method_type | |
69 | { | |
70 | method_euclid = 0, | |
71 | method_binary = 1, | |
72 | method_mixed = 2, | |
73 | }; | |
74 | ||
75 | static const method_type method = | |
76 | boost::has_right_shift_assign<T>::value && boost::has_left_shift_assign<T>::value && boost::has_less<T>::value && boost::has_modulus<T>::value | |
77 | ? method_mixed : | |
78 | boost::has_right_shift_assign<T>::value && boost::has_left_shift_assign<T>::value && boost::has_less<T>::value | |
79 | ? method_binary : method_euclid; | |
80 | }; | |
81 | // | |
82 | // Default gcd_traits just inherits from defaults: | |
83 | // | |
84 | template <class T> | |
85 | struct gcd_traits : public gcd_traits_defaults<T> {}; | |
86 | // | |
87 | // Special handling for polynomials: | |
88 | // | |
89 | namespace tools { | |
90 | template <class T> | |
91 | class polynomial; | |
92 | } | |
93 | ||
94 | template <class T> | |
95 | struct gcd_traits<boost::math::tools::polynomial<T> > : public gcd_traits_defaults<T> | |
96 | { | |
97 | static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; } | |
98 | }; | |
99 | // | |
100 | // Some platforms have fast bitscan operations, that allow us to implement | |
101 | // make_odd much more efficiently: | |
102 | // | |
103 | #if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64)) | |
104 | #pragma intrinsic(_BitScanForward,) | |
105 | template <> | |
106 | struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long> | |
107 | { | |
108 | BOOST_FORCEINLINE static unsigned find_lsb(unsigned long val) | |
109 | { | |
110 | unsigned long result; | |
111 | _BitScanForward(&result, val); | |
112 | return result; | |
113 | } | |
114 | BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val) | |
115 | { | |
116 | unsigned result = find_lsb(val); | |
117 | val >>= result; | |
118 | return result; | |
119 | } | |
120 | }; | |
121 | ||
122 | #ifdef _M_X64 | |
123 | #pragma intrinsic(_BitScanForward64) | |
124 | template <> | |
125 | struct gcd_traits<unsigned __int64> : public gcd_traits_defaults<unsigned __int64> | |
126 | { | |
127 | BOOST_FORCEINLINE static unsigned find_lsb(unsigned __int64 mask) | |
128 | { | |
129 | unsigned long result; | |
130 | _BitScanForward64(&result, mask); | |
131 | return result; | |
132 | } | |
133 | BOOST_FORCEINLINE static unsigned make_odd(unsigned __int64& val) | |
134 | { | |
135 | unsigned result = find_lsb(val); | |
136 | val >>= result; | |
137 | return result; | |
138 | } | |
139 | }; | |
140 | #endif | |
141 | // | |
142 | // Other integer type are trivial adaptations of the above, | |
143 | // this works for signed types too, as by the time these functions | |
144 | // are called, all values are > 0. | |
145 | // | |
146 | template <> struct gcd_traits<long> : public gcd_traits_defaults<long> | |
147 | { BOOST_FORCEINLINE static unsigned make_odd(long& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
148 | template <> struct gcd_traits<unsigned int> : public gcd_traits_defaults<unsigned int> | |
149 | { BOOST_FORCEINLINE static unsigned make_odd(unsigned int& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
150 | template <> struct gcd_traits<int> : public gcd_traits_defaults<int> | |
151 | { BOOST_FORCEINLINE static unsigned make_odd(int& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
152 | template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short> | |
153 | { BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
154 | template <> struct gcd_traits<short> : public gcd_traits_defaults<short> | |
155 | { BOOST_FORCEINLINE static unsigned make_odd(short& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
156 | template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char> | |
157 | { BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
158 | template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char> | |
159 | { BOOST_FORCEINLINE static signed make_odd(signed char& val){ signed result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
160 | template <> struct gcd_traits<char> : public gcd_traits_defaults<char> | |
161 | { BOOST_FORCEINLINE static unsigned make_odd(char& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
162 | template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t> | |
163 | { BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val){ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; | |
164 | #ifdef _M_X64 | |
165 | template <> struct gcd_traits<__int64> : public gcd_traits_defaults<__int64> | |
166 | { BOOST_FORCEINLINE static unsigned make_odd(__int64& val){ unsigned result = gcd_traits<unsigned __int64>::find_lsb(val); val >>= result; return result; } }; | |
167 | #endif | |
168 | ||
169 | #elif defined(BOOST_GCC) || defined(__clang__) || (defined(BOOST_INTEL) && defined(__GNUC__)) | |
170 | ||
171 | template <> | |
172 | struct gcd_traits<unsigned> : public gcd_traits_defaults<unsigned> | |
173 | { | |
174 | BOOST_FORCEINLINE static unsigned find_lsb(unsigned mask) | |
175 | { | |
176 | return __builtin_ctz(mask); | |
177 | } | |
178 | BOOST_FORCEINLINE static unsigned make_odd(unsigned& val) | |
179 | { | |
180 | unsigned result = find_lsb(val); | |
181 | val >>= result; | |
182 | return result; | |
183 | } | |
184 | }; | |
185 | template <> | |
186 | struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long> | |
187 | { | |
188 | BOOST_FORCEINLINE static unsigned find_lsb(unsigned long mask) | |
189 | { | |
190 | return __builtin_ctzl(mask); | |
191 | } | |
192 | BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val) | |
193 | { | |
194 | unsigned result = find_lsb(val); | |
195 | val >>= result; | |
196 | return result; | |
197 | } | |
198 | }; | |
199 | template <> | |
200 | struct gcd_traits<boost::ulong_long_type> : public gcd_traits_defaults<boost::ulong_long_type> | |
201 | { | |
202 | BOOST_FORCEINLINE static unsigned find_lsb(boost::ulong_long_type mask) | |
203 | { | |
204 | return __builtin_ctzll(mask); | |
205 | } | |
206 | BOOST_FORCEINLINE static unsigned make_odd(boost::ulong_long_type& val) | |
207 | { | |
208 | unsigned result = find_lsb(val); | |
209 | val >>= result; | |
210 | return result; | |
211 | } | |
212 | }; | |
213 | // | |
214 | // Other integer type are trivial adaptations of the above, | |
215 | // this works for signed types too, as by the time these functions | |
216 | // are called, all values are > 0. | |
217 | // | |
218 | template <> struct gcd_traits<boost::long_long_type> : public gcd_traits_defaults<boost::long_long_type> | |
219 | { | |
220 | BOOST_FORCEINLINE static unsigned make_odd(boost::long_long_type& val) { unsigned result = gcd_traits<boost::ulong_long_type>::find_lsb(val); val >>= result; return result; } | |
221 | }; | |
222 | template <> struct gcd_traits<long> : public gcd_traits_defaults<long> | |
223 | { | |
224 | BOOST_FORCEINLINE static unsigned make_odd(long& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } | |
225 | }; | |
226 | template <> struct gcd_traits<int> : public gcd_traits_defaults<int> | |
227 | { | |
228 | BOOST_FORCEINLINE static unsigned make_odd(int& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } | |
229 | }; | |
230 | template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short> | |
231 | { | |
232 | BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
233 | }; | |
234 | template <> struct gcd_traits<short> : public gcd_traits_defaults<short> | |
235 | { | |
236 | BOOST_FORCEINLINE static unsigned make_odd(short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
237 | }; | |
238 | template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char> | |
239 | { | |
240 | BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
241 | }; | |
242 | template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char> | |
243 | { | |
244 | BOOST_FORCEINLINE static signed make_odd(signed char& val) { signed result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
245 | }; | |
246 | template <> struct gcd_traits<char> : public gcd_traits_defaults<char> | |
247 | { | |
248 | BOOST_FORCEINLINE static unsigned make_odd(char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
249 | }; | |
250 | template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t> | |
251 | { | |
252 | BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } | |
253 | }; | |
254 | #endif | |
255 | ||
256 | namespace detail | |
257 | { | |
258 | ||
259 | // | |
260 | // The Mixed Binary Euclid Algorithm | |
261 | // Sidi Mohamed Sedjelmaci | |
262 | // Electronic Notes in Discrete Mathematics 35 (2009) 169-176 | |
263 | // | |
264 | template <class T> | |
265 | T mixed_binary_gcd(T u, T v) | |
266 | { | |
267 | using std::swap; | |
268 | if(gcd_traits<T>::less(u, v)) | |
269 | swap(u, v); | |
270 | ||
271 | unsigned shifts = 0; | |
272 | ||
273 | if(!u) | |
274 | return v; | |
275 | if(!v) | |
276 | return u; | |
277 | ||
278 | shifts = (std::min)(gcd_traits<T>::make_odd(u), gcd_traits<T>::make_odd(v)); | |
279 | ||
280 | while(gcd_traits<T>::less(1, v)) | |
281 | { | |
282 | u %= v; | |
283 | v -= u; | |
284 | if(!u) | |
285 | return v << shifts; | |
286 | if(!v) | |
287 | return u << shifts; | |
288 | gcd_traits<T>::make_odd(u); | |
289 | gcd_traits<T>::make_odd(v); | |
290 | if(gcd_traits<T>::less(u, v)) | |
291 | swap(u, v); | |
292 | } | |
293 | return (v == 1 ? v : u) << shifts; | |
294 | } | |
295 | ||
296 | /** Stein gcd (aka 'binary gcd') | |
297 | * | |
298 | * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose | |
299 | */ | |
300 | template <typename SteinDomain> | |
301 | SteinDomain Stein_gcd(SteinDomain m, SteinDomain n) | |
302 | { | |
303 | using std::swap; | |
304 | BOOST_ASSERT(m >= 0); | |
305 | BOOST_ASSERT(n >= 0); | |
306 | if (m == SteinDomain(0)) | |
307 | return n; | |
308 | if (n == SteinDomain(0)) | |
309 | return m; | |
310 | // m > 0 && n > 0 | |
311 | int d_m = gcd_traits<SteinDomain>::make_odd(m); | |
312 | int d_n = gcd_traits<SteinDomain>::make_odd(n); | |
313 | // odd(m) && odd(n) | |
314 | while (m != n) | |
315 | { | |
316 | if (n > m) | |
317 | swap(n, m); | |
318 | m -= n; | |
319 | gcd_traits<SteinDomain>::make_odd(m); | |
320 | } | |
321 | // m == n | |
322 | m <<= (std::min)(d_m, d_n); | |
323 | return m; | |
324 | } | |
325 | ||
326 | ||
327 | /** Euclidean algorithm | |
328 | * | |
329 | * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose | |
330 | * | |
331 | */ | |
332 | template <typename EuclideanDomain> | |
333 | inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) | |
334 | { | |
335 | using std::swap; | |
336 | while (b != EuclideanDomain(0)) | |
337 | { | |
338 | a %= b; | |
339 | swap(a, b); | |
340 | } | |
341 | return a; | |
342 | } | |
343 | ||
344 | ||
345 | template <typename T> | |
346 | inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_mixed, T>::type | |
347 | optimal_gcd_select(T const &a, T const &b) | |
348 | { | |
349 | return detail::mixed_binary_gcd(a, b); | |
350 | } | |
351 | ||
352 | template <typename T> | |
353 | inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_binary, T>::type | |
354 | optimal_gcd_select(T const &a, T const &b) | |
355 | { | |
356 | return detail::Stein_gcd(a, b); | |
357 | } | |
358 | ||
359 | template <typename T> | |
360 | inline BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == gcd_traits<T>::method_euclid, T>::type | |
361 | optimal_gcd_select(T const &a, T const &b) | |
362 | { | |
363 | return detail::Euclid_gcd(a, b); | |
364 | } | |
365 | ||
366 | template <class T> | |
367 | inline T lcm_imp(const T& a, const T& b) | |
368 | { | |
369 | T temp = boost::math::detail::optimal_gcd_select(a, b); | |
370 | #if BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500) | |
371 | return (temp != T(0)) ? T(a / temp * b) : T(0); | |
372 | #else | |
373 | return temp ? T(a / temp * b) : T(0); | |
374 | #endif | |
375 | } | |
376 | ||
377 | } // namespace detail | |
378 | ||
379 | ||
380 | template <typename Integer> | |
381 | inline Integer gcd(Integer const &a, Integer const &b) | |
382 | { | |
383 | return detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b))); | |
384 | } | |
385 | ||
386 | template <typename Integer> | |
387 | inline Integer lcm(Integer const &a, Integer const &b) | |
388 | { | |
389 | return detail::lcm_imp(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b))); | |
390 | } | |
391 | ||
392 | /** | |
393 | * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998 | |
394 | * Chapter 4.5.2, Algorithm C: Greatest common divisor of n integers. | |
395 | * | |
396 | * Knuth counts down from n to zero but we naturally go from first to last. | |
397 | * We also return the termination position because it might be useful to know. | |
398 | * | |
399 | * Partly by quirk, partly by design, this algorithm is defined for n = 1, | |
400 | * because the gcd of {x} is x. It is not defined for n = 0. | |
401 | * | |
402 | * @tparam I Input iterator. | |
403 | * @return The gcd of the range and the iterator position at termination. | |
404 | */ | |
405 | template <typename I> | |
406 | std::pair<typename std::iterator_traits<I>::value_type, I> | |
407 | gcd_range(I first, I last) | |
408 | { | |
409 | BOOST_ASSERT(first != last); | |
410 | typedef typename std::iterator_traits<I>::value_type T; | |
411 | ||
412 | T d = *first++; | |
413 | while (d != T(1) && first != last) | |
414 | { | |
415 | d = gcd(d, *first); | |
416 | first++; | |
417 | } | |
418 | return std::make_pair(d, first); | |
419 | } | |
420 | ||
421 | } // namespace math | |
422 | } // namespace boost | |
423 | ||
424 | #ifdef BOOST_MSVC | |
425 | #pragma warning(pop) | |
426 | #endif | |
427 | ||
428 | #endif // BOOST_MATH_COMMON_FACTOR_RT_HPP |