1 // (C) Copyright Jeremy William Murphy 2016.
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)
7 #ifndef BOOST_MATH_COMMON_FACTOR_RT_HPP
8 #define BOOST_MATH_COMMON_FACTOR_RT_HPP
10 #include <boost/assert.hpp>
11 #include <boost/core/enable_if.hpp>
12 #include <boost/mpl/and.hpp>
13 #include <boost/type_traits.hpp>
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>
23 #if (defined(BOOST_MSVC) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64))
29 #pragma warning(disable:4127 4244) // Conditional expression is constant
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
38 inline static const T& abs(const T& val) { return val; }
41 struct gcd_traits_abs_defaults<T, false>
43 inline static T abs(const T& val)
51 struct gcd_traits_defaults : public gcd_traits_abs_defaults<T>
53 BOOST_FORCEINLINE static unsigned make_odd(T& val)
63 inline static bool less(const T& a, const T& b)
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
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;
82 // Default gcd_traits just inherits from defaults:
85 struct gcd_traits : public gcd_traits_defaults<T> {};
87 // Special handling for polynomials:
95 struct gcd_traits<boost::math::tools::polynomial<T> > : public gcd_traits_defaults<T>
97 static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
100 // Some platforms have fast bitscan operations, that allow us to implement
101 // make_odd much more efficiently:
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,)
106 struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long>
108 BOOST_FORCEINLINE static unsigned find_lsb(unsigned long val)
110 unsigned long result;
111 _BitScanForward(&result, val);
114 BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val)
116 unsigned result = find_lsb(val);
123 #pragma intrinsic(_BitScanForward64)
125 struct gcd_traits<unsigned __int64> : public gcd_traits_defaults<unsigned __int64>
127 BOOST_FORCEINLINE static unsigned find_lsb(unsigned __int64 mask)
129 unsigned long result;
130 _BitScanForward64(&result, mask);
133 BOOST_FORCEINLINE static unsigned make_odd(unsigned __int64& val)
135 unsigned result = find_lsb(val);
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.
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; } };
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; } };
169 #elif defined(BOOST_GCC) || defined(__clang__) || (defined(BOOST_INTEL) && defined(__GNUC__))
172 struct gcd_traits<unsigned> : public gcd_traits_defaults<unsigned>
174 BOOST_FORCEINLINE static unsigned find_lsb(unsigned mask)
176 return __builtin_ctz(mask);
178 BOOST_FORCEINLINE static unsigned make_odd(unsigned& val)
180 unsigned result = find_lsb(val);
186 struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long>
188 BOOST_FORCEINLINE static unsigned find_lsb(unsigned long mask)
190 return __builtin_ctzl(mask);
192 BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val)
194 unsigned result = find_lsb(val);
200 struct gcd_traits<boost::ulong_long_type> : public gcd_traits_defaults<boost::ulong_long_type>
202 BOOST_FORCEINLINE static unsigned find_lsb(boost::ulong_long_type mask)
204 return __builtin_ctzll(mask);
206 BOOST_FORCEINLINE static unsigned make_odd(boost::ulong_long_type& val)
208 unsigned result = find_lsb(val);
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.
218 template <> struct gcd_traits<boost::long_long_type> : public gcd_traits_defaults<boost::long_long_type>
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; }
222 template <> struct gcd_traits<long> : public gcd_traits_defaults<long>
224 BOOST_FORCEINLINE static unsigned make_odd(long& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; }
226 template <> struct gcd_traits<int> : public gcd_traits_defaults<int>
228 BOOST_FORCEINLINE static unsigned make_odd(int& val) { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; }
230 template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short>
232 BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
234 template <> struct gcd_traits<short> : public gcd_traits_defaults<short>
236 BOOST_FORCEINLINE static unsigned make_odd(short& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
238 template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char>
240 BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
242 template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char>
244 BOOST_FORCEINLINE static signed make_odd(signed char& val) { signed result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
246 template <> struct gcd_traits<char> : public gcd_traits_defaults<char>
248 BOOST_FORCEINLINE static unsigned make_odd(char& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
250 template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t>
252 BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val) { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; }
260 // The Mixed Binary Euclid Algorithm
261 // Sidi Mohamed Sedjelmaci
262 // Electronic Notes in Discrete Mathematics 35 (2009) 169-176
265 T mixed_binary_gcd(T u, T v)
268 if(gcd_traits<T>::less(u, v))
278 shifts = (std::min)(gcd_traits<T>::make_odd(u), gcd_traits<T>::make_odd(v));
280 while(gcd_traits<T>::less(1, v))
288 gcd_traits<T>::make_odd(u);
289 gcd_traits<T>::make_odd(v);
290 if(gcd_traits<T>::less(u, v))
293 return (v == 1 ? v : u) << shifts;
296 /** Stein gcd (aka 'binary gcd')
298 * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose
300 template <typename SteinDomain>
301 SteinDomain Stein_gcd(SteinDomain m, SteinDomain n)
304 BOOST_ASSERT(m >= 0);
305 BOOST_ASSERT(n >= 0);
306 if (m == SteinDomain(0))
308 if (n == SteinDomain(0))
311 int d_m = gcd_traits<SteinDomain>::make_odd(m);
312 int d_n = gcd_traits<SteinDomain>::make_odd(n);
319 gcd_traits<SteinDomain>::make_odd(m);
322 m <<= (std::min)(d_m, d_n);
327 /** Euclidean algorithm
329 * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose
332 template <typename EuclideanDomain>
333 inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b)
336 while (b != EuclideanDomain(0))
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)
349 return detail::mixed_binary_gcd(a, b);
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)
356 return detail::Stein_gcd(a, b);
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)
363 return detail::Euclid_gcd(a, b);
367 inline T lcm_imp(const T& a, const T& b)
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);
373 return temp ? T(a / temp * b) : T(0);
377 } // namespace detail
380 template <typename Integer>
381 inline Integer gcd(Integer const &a, Integer const &b)
383 return detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
386 template <typename Integer>
387 inline Integer lcm(Integer const &a, Integer const &b)
389 return detail::lcm_imp(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
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.
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.
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.
402 * @tparam I Input iterator.
403 * @return The gcd of the range and the iterator position at termination.
405 template <typename I>
406 std::pair<typename std::iterator_traits<I>::value_type, I>
407 gcd_range(I first, I last)
409 BOOST_ASSERT(first != last);
410 typedef typename std::iterator_traits<I>::value_type T;
413 while (d != T(1) && first != last)
418 return std::make_pair(d, first);
428 #endif // BOOST_MATH_COMMON_FACTOR_RT_HPP