diff options
Diffstat (limited to 'boost/integer/common_factor_rt.hpp')
| -rw-r--r-- | boost/integer/common_factor_rt.hpp | 898 |
1 files changed, 508 insertions, 390 deletions
diff --git a/boost/integer/common_factor_rt.hpp b/boost/integer/common_factor_rt.hpp index c2b54db8e8..341b316501 100644 --- a/boost/integer/common_factor_rt.hpp +++ b/boost/integer/common_factor_rt.hpp @@ -1,454 +1,572 @@ -// Boost common_factor_rt.hpp header file ----------------------------------// +// (C) Copyright Jeremy William Murphy 2016. -// (C) Copyright Daryle Walker and Paul Moore 2001-2002. Permission to copy, -// use, modify, sell and distribute this software is granted provided this -// copyright notice appears in all copies. This software is provided "as is" -// without express or implied warranty, and with no claim as to its suitability -// for any purpose. - -// boostinspect:nolicense (don't complain about the lack of a Boost license) -// (Paul Moore hasn't been in contact for years, so there's no way to change the -// license.) - -// See http://www.boost.org for updates, documentation, and revision history. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_INTEGER_COMMON_FACTOR_RT_HPP #define BOOST_INTEGER_COMMON_FACTOR_RT_HPP -#include <boost/integer_fwd.hpp> // self include +#include <boost/assert.hpp> +#include <boost/core/enable_if.hpp> #include <boost/config.hpp> // for BOOST_NESTED_TEMPLATE, etc. #include <boost/limits.hpp> // for std::numeric_limits #include <climits> // for CHAR_MIN #include <boost/detail/workaround.hpp> +#include <iterator> +#include <algorithm> +#include <limits> +#ifndef BOOST_NO_CXX11_HDR_TYPE_TRAITS +#include <type_traits> +#endif +#ifdef BOOST_NO_CXX11_HDR_FUNCTIONAL +#include <functional> +#endif + +#if ((defined(BOOST_MSVC) && (BOOST_MSVC >= 1600)) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64)) +#include <intrin.h> +#endif #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable:4127 4244) // Conditional expression is constant #endif -namespace boost -{ -namespace integer -{ - - -// Forward declarations for function templates -----------------------------// - -template < typename IntegerType > - IntegerType gcd( IntegerType const &a, IntegerType const &b ); - -template < typename IntegerType > - IntegerType lcm( IntegerType const &a, IntegerType const &b ); - - -// Greatest common divisor evaluator class declaration ---------------------// - -template < typename IntegerType > -class gcd_evaluator -{ -public: - // Types - typedef IntegerType result_type, first_argument_type, second_argument_type; - - // Function object interface - result_type operator ()( first_argument_type const &a, - second_argument_type const &b ) const; - -}; // boost::integer::gcd_evaluator - - -// Least common multiple evaluator class declaration -----------------------// - -template < typename IntegerType > -class lcm_evaluator -{ -public: - // Types - typedef IntegerType result_type, first_argument_type, second_argument_type; - - // Function object interface - result_type operator ()( first_argument_type const &a, - second_argument_type const &b ) const; - -}; // boost::integer::lcm_evaluator - - -// Implementation details --------------------------------------------------// - -namespace detail -{ - // Greatest common divisor for rings (including unsigned integers) - template < typename RingType > - RingType - gcd_euclidean - ( - RingType a, - RingType b - ) - { - // Avoid repeated construction - #ifndef __BORLANDC__ - RingType const zero = static_cast<RingType>( 0 ); - #else - RingType zero = static_cast<RingType>( 0 ); - #endif - - // Reduce by GCD-remainder property [GCD(a,b) == GCD(b,a MOD b)] - while ( true ) - { - if ( a == zero ) - return b; - b %= a; - - if ( b == zero ) - return a; - a %= b; - } - } - - // Greatest common divisor for (signed) integers - template < typename IntegerType > - inline - IntegerType - gcd_integer - ( - IntegerType const & a, - IntegerType const & b - ) - { - // Avoid repeated construction - IntegerType const zero = static_cast<IntegerType>( 0 ); - IntegerType const result = gcd_euclidean( a, b ); - - return ( result < zero ) ? static_cast<IntegerType>(-result) : result; - } - - // Greatest common divisor for unsigned binary integers - template < typename BuiltInUnsigned > - BuiltInUnsigned - gcd_binary - ( - BuiltInUnsigned u, - BuiltInUnsigned v - ) - { - if ( u && v ) - { - // Shift out common factors of 2 - unsigned shifts = 0; - - while ( !(u & 1u) && !(v & 1u) ) - { - ++shifts; - u >>= 1; - v >>= 1; - } +#if !defined(BOOST_NO_CXX11_HDR_TYPE_TRAITS) && !defined(BOOST_NO_CXX11_NOEXCEPT) +#define BOOST_GCD_NOEXCEPT(T) noexcept(std::is_arithmetic<T>::value) +#else +#define BOOST_GCD_NOEXCEPT(T) +#endif - // Start with the still-even one, if any - BuiltInUnsigned r[] = { u, v }; - unsigned which = static_cast<bool>( u & 1u ); +namespace boost { + + template <class I> + class rational; + + namespace integer { + + namespace gcd_detail{ + + // + // some helper functions which really should be constexpr already, but sadly aren't: + // +#ifndef BOOST_NO_CXX14_CONSTEXPR + template <class T> + inline constexpr T constexpr_min(T const& a, T const& b) BOOST_GCD_NOEXCEPT(T) + { + return a < b ? a : b; + } + template <class T> + inline constexpr auto constexpr_swap(T&a, T& b) BOOST_GCD_NOEXCEPT(T) -> decltype(a.swap(b)) + { + return a.swap(b); + } + template <class T, class U> + inline constexpr void constexpr_swap(T&a, U& b...) BOOST_GCD_NOEXCEPT(T) + { + T t(static_cast<T&&>(a)); + a = static_cast<T&&>(b); + b = static_cast<T&&>(t); + } +#else + template <class T> + inline T constexpr_min(T const& a, T const& b) BOOST_GCD_NOEXCEPT(T) + { + return a < b ? a : b; + } + template <class T> + inline void constexpr_swap(T&a, T& b) BOOST_GCD_NOEXCEPT(T) + { + using std::swap; + swap(a, b); + } +#endif - // Whittle down the values via their differences - do + template <class T, bool a = +#ifndef BOOST_NO_CXX11_HDR_TYPE_TRAITS + std::is_unsigned<T>::value || +#endif + (std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_signed)> + struct gcd_traits_abs_defaults + { + inline static BOOST_CXX14_CONSTEXPR const T& abs(const T& val) BOOST_GCD_NOEXCEPT(T) { return val; } + }; + template <class T> + struct gcd_traits_abs_defaults<T, false> + { + inline static T BOOST_CXX14_CONSTEXPR abs(const T& val) BOOST_GCD_NOEXCEPT(T) + { + // This sucks, but std::abs is not constexpr :( + return val < T(0) ? -val : val; + } + }; + + enum method_type + { + method_euclid = 0, + method_binary = 1, + method_mixed = 2 + }; + + struct any_convert + { + template <class T> + any_convert(const T&); + }; + + struct unlikely_size + { + char buf[9973]; + }; + + unlikely_size operator <<= (any_convert, any_convert); + unlikely_size operator >>= (any_convert, any_convert); + + template <class T> + struct gcd_traits_defaults : public gcd_traits_abs_defaults<T> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(T& val) BOOST_GCD_NOEXCEPT(T) + { + unsigned r = 0; + while(0 == (val & 1u)) { -#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) - while ( !(r[ which ] & 1u) ) - { - r[ which ] = (r[which] >> 1); - } +#ifdef _MSC_VER // VC++ can't handle operator >>= in constexpr code for some reason + val = val >> 1; #else - // Remove factors of two from the even one - while ( !(r[ which ] & 1u) ) - { - r[ which ] >>= 1; - } + val >>= 1; #endif - - // Replace the larger of the two with their difference - if ( r[!which] > r[which] ) - { - which ^= 1u; - } - - r[ which ] -= r[ !which ]; + ++r; } - while ( r[which] ); - - // Shift-in the common factor of 2 to the residues' GCD - return r[ !which ] << shifts; - } - else - { - // At least one input is zero, return the other - // (adding since zero is the additive identity) - // or zero if both are zero. - return u + v; - } - } - - // Least common multiple for rings (including unsigned integers) - template < typename RingType > - inline - RingType - lcm_euclidean - ( - RingType const & a, - RingType const & b - ) - { - RingType const zero = static_cast<RingType>( 0 ); - RingType const temp = gcd_euclidean( a, b ); - - return ( temp != zero ) ? ( a / temp * b ) : zero; - } - - // Least common multiple for (signed) integers - template < typename IntegerType > - inline - IntegerType - lcm_integer - ( - IntegerType const & a, - IntegerType const & b - ) - { - // Avoid repeated construction - IntegerType const zero = static_cast<IntegerType>( 0 ); - IntegerType const result = lcm_euclidean( a, b ); - - return ( result < zero ) ? static_cast<IntegerType>(-result) : result; - } - - // Function objects to find the best way of computing GCD or LCM -#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - template < typename T, bool IsSpecialized, bool IsSigned > - struct gcd_optimal_evaluator_helper_t - { - T operator ()( T const &a, T const &b ) - { - return gcd_euclidean( a, b ); - } - }; - - template < typename T > - struct gcd_optimal_evaluator_helper_t< T, true, true > - { - T operator ()( T const &a, T const &b ) - { - return gcd_integer( a, b ); - } - }; - - template < typename T > - struct gcd_optimal_evaluator - { - T operator ()( T const &a, T const &b ) - { - typedef ::std::numeric_limits<T> limits_type; - - typedef gcd_optimal_evaluator_helper_t<T, - limits_type::is_specialized, limits_type::is_signed> helper_type; - - helper_type solver; - - return solver( a, b ); - } - }; -#else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - template < typename T > - struct gcd_optimal_evaluator - { - T operator ()( T const &a, T const &b ) - { - return gcd_integer( a, b ); - } - }; + return r; + } + inline static BOOST_CXX14_CONSTEXPR bool less(const T& a, const T& b) BOOST_GCD_NOEXCEPT(T) + { + return a < b; + } + + static T& get_value(); + +#ifndef BOOST_NO_SFINAE + static const bool has_operator_left_shift_equal = sizeof(get_value() <<= 2) != sizeof(unlikely_size); + static const bool has_operator_right_shift_equal = sizeof(get_value() >>= 2) != sizeof(unlikely_size); +#else + static const bool has_operator_left_shift_equal = true; + static const bool has_operator_right_shift_equal = true; #endif - - // Specialize for the built-in integers -#define BOOST_PRIVATE_GCD_UF( Ut ) \ - template < > struct gcd_optimal_evaluator<Ut> \ - { Ut operator ()( Ut a, Ut b ) const { return gcd_binary( a, b ); } } - - BOOST_PRIVATE_GCD_UF( unsigned char ); - BOOST_PRIVATE_GCD_UF( unsigned short ); - BOOST_PRIVATE_GCD_UF( unsigned ); - BOOST_PRIVATE_GCD_UF( unsigned long ); - -#ifdef BOOST_HAS_LONG_LONG - BOOST_PRIVATE_GCD_UF( boost::ulong_long_type ); -#elif defined(BOOST_HAS_MS_INT64) - BOOST_PRIVATE_GCD_UF( unsigned __int64 ); + static const method_type method = std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::is_integer && has_operator_left_shift_equal && has_operator_right_shift_equal ? method_mixed : method_euclid; + }; + // + // Default gcd_traits just inherits from defaults: + // + template <class T> + struct gcd_traits : public gcd_traits_defaults<T> {}; + + // + // Some platforms have fast bitscan operations, that allow us to implement + // make_odd much more efficiently, unfortunately we can't use these if we want + // the functions to be constexpr as the compiler intrinsics aren't constexpr. + // +#if defined(BOOST_NO_CXX14_CONSTEXPR) && ((defined(BOOST_MSVC) && (BOOST_MSVC >= 1600)) || (defined(__clang__) && defined(__c2__)) || (defined(BOOST_INTEL) && defined(_MSC_VER))) && (defined(_M_IX86) || defined(_M_X64)) +#pragma intrinsic(_BitScanForward,) + template <> + struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long> + { + BOOST_FORCEINLINE static unsigned find_lsb(unsigned long val) BOOST_NOEXCEPT + { + unsigned long result; + _BitScanForward(&result, val); + return result; + } + BOOST_FORCEINLINE static unsigned make_odd(unsigned long& val) BOOST_NOEXCEPT + { + unsigned result = find_lsb(val); + val >>= result; + return result; + } + }; + +#ifdef _M_X64 +#pragma intrinsic(_BitScanForward64) + template <> + struct gcd_traits<unsigned __int64> : public gcd_traits_defaults<unsigned __int64> + { + BOOST_FORCEINLINE static unsigned find_lsb(unsigned __int64 mask) BOOST_NOEXCEPT + { + unsigned long result; + _BitScanForward64(&result, mask); + return result; + } + BOOST_FORCEINLINE static unsigned make_odd(unsigned __int64& val) BOOST_NOEXCEPT + { + unsigned result = find_lsb(val); + val >>= result; + return result; + } + }; #endif - -#if CHAR_MIN == 0 - BOOST_PRIVATE_GCD_UF( char ); // char is unsigned + // + // Other integer type are trivial adaptations of the above, + // this works for signed types too, as by the time these functions + // are called, all values are > 0. + // + template <> struct gcd_traits<long> : public gcd_traits_defaults<long> + { BOOST_FORCEINLINE static unsigned make_odd(long& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<unsigned int> : public gcd_traits_defaults<unsigned int> + { BOOST_FORCEINLINE static unsigned make_odd(unsigned int& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<int> : public gcd_traits_defaults<int> + { BOOST_FORCEINLINE static unsigned make_odd(int& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short> + { BOOST_FORCEINLINE static unsigned make_odd(unsigned short& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<short> : public gcd_traits_defaults<short> + { BOOST_FORCEINLINE static unsigned make_odd(short& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char> + { BOOST_FORCEINLINE static unsigned make_odd(unsigned char& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char> + { BOOST_FORCEINLINE static signed make_odd(signed char& val)BOOST_NOEXCEPT{ signed result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; + template <> struct gcd_traits<char> : public gcd_traits_defaults<char> + { BOOST_FORCEINLINE static unsigned make_odd(char& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; +#ifndef BOOST_NO_INTRINSIC_WCHAR_T + template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t> + { BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } }; #endif - -#undef BOOST_PRIVATE_GCD_UF - -#define BOOST_PRIVATE_GCD_SF( St, Ut ) \ - template < > struct gcd_optimal_evaluator<St> \ - { St operator ()( St a, St b ) const { Ut const a_abs = \ - static_cast<Ut>( a < 0 ? -a : +a ), b_abs = static_cast<Ut>( \ - b < 0 ? -b : +b ); return static_cast<St>( \ - gcd_optimal_evaluator<Ut>()(a_abs, b_abs) ); } } - - BOOST_PRIVATE_GCD_SF( signed char, unsigned char ); - BOOST_PRIVATE_GCD_SF( short, unsigned short ); - BOOST_PRIVATE_GCD_SF( int, unsigned ); - BOOST_PRIVATE_GCD_SF( long, unsigned long ); - -#if CHAR_MIN < 0 - BOOST_PRIVATE_GCD_SF( char, unsigned char ); // char is signed +#ifdef _M_X64 + template <> struct gcd_traits<__int64> : public gcd_traits_defaults<__int64> + { BOOST_FORCEINLINE static unsigned make_odd(__int64& val)BOOST_NOEXCEPT{ unsigned result = gcd_traits<unsigned __int64>::find_lsb(val); val >>= result; return result; } }; #endif -#ifdef BOOST_HAS_LONG_LONG - BOOST_PRIVATE_GCD_SF( boost::long_long_type, boost::ulong_long_type ); -#elif defined(BOOST_HAS_MS_INT64) - BOOST_PRIVATE_GCD_SF( __int64, unsigned __int64 ); +#elif defined(BOOST_GCC) || defined(__clang__) || (defined(BOOST_INTEL) && defined(__GNUC__)) + + template <> + struct gcd_traits<unsigned> : public gcd_traits_defaults<unsigned> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned find_lsb(unsigned mask)BOOST_NOEXCEPT + { + return __builtin_ctz(mask); + } + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(unsigned& val)BOOST_NOEXCEPT + { + unsigned result = find_lsb(val); + val >>= result; + return result; + } + }; + template <> + struct gcd_traits<unsigned long> : public gcd_traits_defaults<unsigned long> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned find_lsb(unsigned long mask)BOOST_NOEXCEPT + { + return __builtin_ctzl(mask); + } + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(unsigned long& val)BOOST_NOEXCEPT + { + unsigned result = find_lsb(val); + val >>= result; + return result; + } + }; + template <> + struct gcd_traits<boost::ulong_long_type> : public gcd_traits_defaults<boost::ulong_long_type> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned find_lsb(boost::ulong_long_type mask)BOOST_NOEXCEPT + { + return __builtin_ctzll(mask); + } + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(boost::ulong_long_type& val)BOOST_NOEXCEPT + { + unsigned result = find_lsb(val); + val >>= result; + return result; + } + }; + // + // Other integer type are trivial adaptations of the above, + // this works for signed types too, as by the time these functions + // are called, all values are > 0. + // + template <> struct gcd_traits<boost::long_long_type> : public gcd_traits_defaults<boost::long_long_type> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(boost::long_long_type& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<boost::ulong_long_type>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<long> : public gcd_traits_defaults<long> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(long& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<int> : public gcd_traits_defaults<int> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(int& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned long>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<unsigned short> : public gcd_traits_defaults<unsigned short> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(unsigned short& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<short> : public gcd_traits_defaults<short> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(short& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<unsigned char> : public gcd_traits_defaults<unsigned char> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(unsigned char& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<signed char> : public gcd_traits_defaults<signed char> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR signed make_odd(signed char& val)BOOST_NOEXCEPT { signed result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; + template <> struct gcd_traits<char> : public gcd_traits_defaults<char> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(char& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; +#ifndef BOOST_NO_INTRINSIC_WCHAR_T + template <> struct gcd_traits<wchar_t> : public gcd_traits_defaults<wchar_t> + { + BOOST_FORCEINLINE static BOOST_CXX14_CONSTEXPR unsigned make_odd(wchar_t& val)BOOST_NOEXCEPT { unsigned result = gcd_traits<unsigned>::find_lsb(val); val >>= result; return result; } + }; #endif - -#undef BOOST_PRIVATE_GCD_SF - -#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - template < typename T, bool IsSpecialized, bool IsSigned > - struct lcm_optimal_evaluator_helper_t +#endif + // + // The Mixed Binary Euclid Algorithm + // Sidi Mohamed Sedjelmaci + // Electronic Notes in Discrete Mathematics 35 (2009) 169-176 + // + template <class T> + BOOST_CXX14_CONSTEXPR T mixed_binary_gcd(T u, T v) BOOST_GCD_NOEXCEPT(T) + { + if(gcd_traits<T>::less(u, v)) + constexpr_swap(u, v); + + unsigned shifts = 0; + + if(u == T(0)) + return v; + if(v == T(0)) + return u; + + shifts = constexpr_min(gcd_traits<T>::make_odd(u), gcd_traits<T>::make_odd(v)); + + while(gcd_traits<T>::less(1, v)) + { + u %= v; + v -= u; + if(u == T(0)) + return v << shifts; + if(v == T(0)) + return u << shifts; + gcd_traits<T>::make_odd(u); + gcd_traits<T>::make_odd(v); + if(gcd_traits<T>::less(u, v)) + constexpr_swap(u, v); + } + return (v == 1 ? v : u) << shifts; + } + + /** Stein gcd (aka 'binary gcd') + * + * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose + */ + template <typename SteinDomain> + BOOST_CXX14_CONSTEXPR SteinDomain Stein_gcd(SteinDomain m, SteinDomain n) BOOST_GCD_NOEXCEPT(SteinDomain) { - T operator ()( T const &a, T const &b ) + BOOST_ASSERT(m >= 0); + BOOST_ASSERT(n >= 0); + if (m == SteinDomain(0)) + return n; + if (n == SteinDomain(0)) + return m; + // m > 0 && n > 0 + int d_m = gcd_traits<SteinDomain>::make_odd(m); + int d_n = gcd_traits<SteinDomain>::make_odd(n); + // odd(m) && odd(n) + while (m != n) { - return lcm_euclidean( a, b ); + if (n > m) + constexpr_swap(n, m); + m -= n; + gcd_traits<SteinDomain>::make_odd(m); } - }; + // m == n + m <<= constexpr_min(d_m, d_n); + return m; + } - template < typename T > - struct lcm_optimal_evaluator_helper_t< T, true, true > + + /** Euclidean algorithm + * + * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose + * + */ + template <typename EuclideanDomain> + inline BOOST_CXX14_CONSTEXPR EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) BOOST_GCD_NOEXCEPT(EuclideanDomain) { - T operator ()( T const &a, T const &b ) + while (b != EuclideanDomain(0)) { - return lcm_integer( a, b ); + a %= b; + constexpr_swap(a, b); } - }; - - template < typename T > - struct lcm_optimal_evaluator - { - T operator ()( T const &a, T const &b ) - { - typedef ::std::numeric_limits<T> limits_type; - - typedef lcm_optimal_evaluator_helper_t<T, - limits_type::is_specialized, limits_type::is_signed> helper_type; + return a; + } - helper_type solver; - return solver( a, b ); - } - }; -#else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS - template < typename T > - struct lcm_optimal_evaluator + template <typename T> + inline BOOST_CXX14_CONSTEXPR BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == method_mixed, T>::type + optimal_gcd_select(T const &a, T const &b) BOOST_GCD_NOEXCEPT(T) { - T operator ()( T const &a, T const &b ) - { - return lcm_integer( a, b ); - } - }; -#endif + return gcd_detail::mixed_binary_gcd(a, b); + } - // Functions to find the GCD or LCM in the best way - template < typename T > - inline - T - gcd_optimal - ( - T const & a, - T const & b - ) + template <typename T> + inline BOOST_CXX14_CONSTEXPR BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == method_binary, T>::type + optimal_gcd_select(T const &a, T const &b) BOOST_GCD_NOEXCEPT(T) { - gcd_optimal_evaluator<T> solver; - - return solver( a, b ); + return gcd_detail::Stein_gcd(a, b); } - template < typename T > - inline - T - lcm_optimal - ( - T const & a, - T const & b - ) + template <typename T> + inline BOOST_CXX14_CONSTEXPR BOOST_DEDUCED_TYPENAME enable_if_c<gcd_traits<T>::method == method_euclid, T>::type + optimal_gcd_select(T const &a, T const &b) BOOST_GCD_NOEXCEPT(T) { - lcm_optimal_evaluator<T> solver; - - return solver( a, b ); + return gcd_detail::Euclid_gcd(a, b); } -} // namespace detail + template <class T> + inline BOOST_CXX14_CONSTEXPR T lcm_imp(const T& a, const T& b) BOOST_GCD_NOEXCEPT(T) + { + T temp = boost::integer::gcd_detail::optimal_gcd_select(a, b); +#if BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500) + return (temp != T(0)) ? T(a / temp * b) : T(0); +#else + return temp != T(0) ? T(a / temp * b) : T(0); +#endif + } +} // namespace detail -// Greatest common divisor evaluator member function definition ------------// -template < typename IntegerType > -inline -typename gcd_evaluator<IntegerType>::result_type -gcd_evaluator<IntegerType>::operator () -( - first_argument_type const & a, - second_argument_type const & b -) const +template <typename Integer> +inline BOOST_CXX14_CONSTEXPR Integer gcd(Integer const &a, Integer const &b) BOOST_GCD_NOEXCEPT(Integer) { - return detail::gcd_optimal( a, b ); + if(a == (std::numeric_limits<Integer>::min)()) + return a == static_cast<Integer>(0) ? gcd_detail::gcd_traits<Integer>::abs(b) : boost::integer::gcd(static_cast<Integer>(a % b), b); + else if (b == (std::numeric_limits<Integer>::min)()) + return b == static_cast<Integer>(0) ? gcd_detail::gcd_traits<Integer>::abs(a) : boost::integer::gcd(a, static_cast<Integer>(b % a)); + return gcd_detail::optimal_gcd_select(static_cast<Integer>(gcd_detail::gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_detail::gcd_traits<Integer>::abs(b))); } - -// Least common multiple evaluator member function definition --------------// - -template < typename IntegerType > -inline -typename lcm_evaluator<IntegerType>::result_type -lcm_evaluator<IntegerType>::operator () -( - first_argument_type const & a, - second_argument_type const & b -) const +template <typename Integer> +inline BOOST_CXX14_CONSTEXPR Integer lcm(Integer const &a, Integer const &b) BOOST_GCD_NOEXCEPT(Integer) { - return detail::lcm_optimal( a, b ); + return gcd_detail::lcm_imp(static_cast<Integer>(gcd_detail::gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_detail::gcd_traits<Integer>::abs(b))); +} +#ifndef BOOST_NO_CXX11_VARIADIC_TEMPLATES +// +// This looks slightly odd, but the variadic forms must have 3 or more arguments, and the variadic argument pack may be empty. +// This matters not at all for most compilers, but Oracle C++ selects the wrong overload in the 2-arg case unless we do this. +// +template <typename Integer, typename... Args> +inline BOOST_CXX14_CONSTEXPR Integer gcd(Integer const &a, Integer const &b, const Integer& c, Args const&... args) BOOST_GCD_NOEXCEPT(Integer) +{ + Integer t = gcd(b, c, args...); + return t == 1 ? 1 : gcd(a, t); } - -// Greatest common divisor and least common multiple function definitions --// - -template < typename IntegerType > -inline -IntegerType -gcd -( - IntegerType const & a, - IntegerType const & b -) +template <typename Integer, typename... Args> +inline BOOST_CXX14_CONSTEXPR Integer lcm(Integer const &a, Integer const &b, Integer const& c, Args const&... args) BOOST_GCD_NOEXCEPT(Integer) { - gcd_evaluator<IntegerType> solver; + return lcm(a, lcm(b, c, args...)); +} +#endif +// +// Special handling for rationals: +// +template <typename Integer> +inline typename boost::enable_if_c<std::numeric_limits<Integer>::is_specialized, boost::rational<Integer> >::type gcd(boost::rational<Integer> const &a, boost::rational<Integer> const &b) +{ + return boost::rational<Integer>(static_cast<Integer>(gcd(a.numerator(), b.numerator())), static_cast<Integer>(lcm(a.denominator(), b.denominator()))); +} - return solver( a, b ); +template <typename Integer> +inline typename boost::enable_if_c<std::numeric_limits<Integer>::is_specialized, boost::rational<Integer> >::type lcm(boost::rational<Integer> const &a, boost::rational<Integer> const &b) +{ + return boost::rational<Integer>(static_cast<Integer>(lcm(a.numerator(), b.numerator())), static_cast<Integer>(gcd(a.denominator(), b.denominator()))); +} +/** + * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998 + * Chapter 4.5.2, Algorithm C: Greatest common divisor of n integers. + * + * Knuth counts down from n to zero but we naturally go from first to last. + * We also return the termination position because it might be useful to know. + * + * Partly by quirk, partly by design, this algorithm is defined for n = 1, + * because the gcd of {x} is x. It is not defined for n = 0. + * + * @tparam I Input iterator. + * @return The gcd of the range and the iterator position at termination. + */ +template <typename I> +std::pair<typename std::iterator_traits<I>::value_type, I> +gcd_range(I first, I last) BOOST_GCD_NOEXCEPT(I) +{ + BOOST_ASSERT(first != last); + typedef typename std::iterator_traits<I>::value_type T; + + T d = *first++; + while (d != T(1) && first != last) + { + d = gcd(d, *first); + first++; + } + return std::make_pair(d, first); +} +template <typename I> +std::pair<typename std::iterator_traits<I>::value_type, I> +lcm_range(I first, I last) BOOST_GCD_NOEXCEPT(I) +{ + BOOST_ASSERT(first != last); + typedef typename std::iterator_traits<I>::value_type T; + + T d = *first++; + while (d != T(1) && first != last) + { + d = lcm(d, *first); + first++; + } + return std::make_pair(d, first); } template < typename IntegerType > -inline -IntegerType -lcm -( - IntegerType const & a, - IntegerType const & b -) +class gcd_evaluator +#ifdef BOOST_NO_CXX11_HDR_FUNCTIONAL + : public std::binary_function<IntegerType, IntegerType, IntegerType> +#endif { - lcm_evaluator<IntegerType> solver; - - return solver( a, b ); -} +public: +#ifndef BOOST_NO_CXX11_HDR_FUNCTIONAL + typedef IntegerType first_argument_type; + typedef IntegerType second_argument_type; + typedef IntegerType result_type; +#endif + IntegerType operator()(IntegerType const &a, IntegerType const &b)const + { + return boost::integer::gcd(a, b); + } +}; +template < typename IntegerType > +class lcm_evaluator +#ifdef BOOST_NO_CXX11_HDR_FUNCTIONAL + : public std::binary_function<IntegerType, IntegerType, IntegerType> +#endif +{ +public: +#ifndef BOOST_NO_CXX11_HDR_FUNCTIONAL + typedef IntegerType first_argument_type; + typedef IntegerType second_argument_type; + typedef IntegerType result_type; +#endif + IntegerType operator()(IntegerType const &a, IntegerType const &b)const + { + return boost::integer::lcm(a, b); + } +}; } // namespace integer } // namespace boost |