Chris@16: //  Boost common_factor_rt.hpp header file  ----------------------------------//
Chris@16: 
Chris@16: //  (C) Copyright Daryle Walker and Paul Moore 2001-2002.  Permission to copy,
Chris@16: //  use, modify, sell and distribute this software is granted provided this
Chris@16: //  copyright notice appears in all copies.  This software is provided "as is"
Chris@16: //  without express or implied warranty, and with no claim as to its suitability
Chris@16: //  for any purpose. 
Chris@16: 
Chris@16: // boostinspect:nolicense (don't complain about the lack of a Boost license)
Chris@16: // (Paul Moore hasn't been in contact for years, so there's no way to change the
Chris@16: // license.)
Chris@16: 
Chris@16: //  See http://www.boost.org for updates, documentation, and revision history. 
Chris@16: 
Chris@16: #ifndef BOOST_MATH_COMMON_FACTOR_RT_HPP
Chris@16: #define BOOST_MATH_COMMON_FACTOR_RT_HPP
Chris@16: 
Chris@16: #include <boost/math_fwd.hpp>  // self include
Chris@16: 
Chris@16: #include <boost/config.hpp>  // for BOOST_NESTED_TEMPLATE, etc.
Chris@16: #include <boost/limits.hpp>  // for std::numeric_limits
Chris@16: #include <climits>           // for CHAR_MIN
Chris@16: #include <boost/detail/workaround.hpp>
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push)
Chris@16: #pragma warning(disable:4127 4244)  // Conditional expression is constant
Chris@16: #endif
Chris@16: 
Chris@16: namespace boost
Chris@16: {
Chris@16: namespace math
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: //  Forward declarations for function templates  -----------------------------//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16:     IntegerType  gcd( IntegerType const &a, IntegerType const &b );
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16:     IntegerType  lcm( IntegerType const &a, IntegerType const &b );
Chris@16: 
Chris@16: 
Chris@16: //  Greatest common divisor evaluator class declaration  ---------------------//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: class gcd_evaluator
Chris@16: {
Chris@16: public:
Chris@16:     // Types
Chris@16:     typedef IntegerType  result_type, first_argument_type, second_argument_type;
Chris@16: 
Chris@16:     // Function object interface
Chris@16:     result_type  operator ()( first_argument_type const &a,
Chris@16:      second_argument_type const &b ) const;
Chris@16: 
Chris@16: };  // boost::math::gcd_evaluator
Chris@16: 
Chris@16: 
Chris@16: //  Least common multiple evaluator class declaration  -----------------------//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: class lcm_evaluator
Chris@16: {
Chris@16: public:
Chris@16:     // Types
Chris@16:     typedef IntegerType  result_type, first_argument_type, second_argument_type;
Chris@16: 
Chris@16:     // Function object interface
Chris@16:     result_type  operator ()( first_argument_type const &a,
Chris@16:      second_argument_type const &b ) const;
Chris@16: 
Chris@16: };  // boost::math::lcm_evaluator
Chris@16: 
Chris@16: 
Chris@16: //  Implementation details  --------------------------------------------------//
Chris@16: 
Chris@16: namespace detail
Chris@16: {
Chris@16:     // Greatest common divisor for rings (including unsigned integers)
Chris@16:     template < typename RingType >
Chris@16:     RingType
Chris@16:     gcd_euclidean
Chris@16:     (
Chris@16:         RingType a,
Chris@16:         RingType b
Chris@16:     )
Chris@16:     {
Chris@16:         // Avoid repeated construction
Chris@16:         #ifndef __BORLANDC__
Chris@16:         RingType const  zero = static_cast<RingType>( 0 );
Chris@16:         #else
Chris@16:         RingType  zero = static_cast<RingType>( 0 );
Chris@16:         #endif
Chris@16: 
Chris@16:         // Reduce by GCD-remainder property [GCD(a,b) == GCD(b,a MOD b)]
Chris@16:         while ( true )
Chris@16:         {
Chris@16:             if ( a == zero )
Chris@16:                 return b;
Chris@16:             b %= a;
Chris@16: 
Chris@16:             if ( b == zero )
Chris@16:                 return a;
Chris@16:             a %= b;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     // Greatest common divisor for (signed) integers
Chris@16:     template < typename IntegerType >
Chris@16:     inline
Chris@16:     IntegerType
Chris@16:     gcd_integer
Chris@16:     (
Chris@16:         IntegerType const &  a,
Chris@16:         IntegerType const &  b
Chris@16:     )
Chris@16:     {
Chris@16:         // Avoid repeated construction
Chris@16:         IntegerType const  zero = static_cast<IntegerType>( 0 );
Chris@16:         IntegerType const  result = gcd_euclidean( a, b );
Chris@16: 
Chris@16:         return ( result < zero ) ? static_cast<IntegerType>(-result) : result;
Chris@16:     }
Chris@16: 
Chris@16:     // Greatest common divisor for unsigned binary integers
Chris@16:     template < typename BuiltInUnsigned >
Chris@16:     BuiltInUnsigned
Chris@16:     gcd_binary
Chris@16:     (
Chris@16:         BuiltInUnsigned  u,
Chris@16:         BuiltInUnsigned  v
Chris@16:     )
Chris@16:     {
Chris@16:         if ( u && v )
Chris@16:         {
Chris@16:             // Shift out common factors of 2
Chris@16:             unsigned  shifts = 0;
Chris@16: 
Chris@16:             while ( !(u & 1u) && !(v & 1u) )
Chris@16:             {
Chris@16:                 ++shifts;
Chris@16:                 u >>= 1;
Chris@16:                 v >>= 1;
Chris@16:             }
Chris@16: 
Chris@16:             // Start with the still-even one, if any
Chris@16:             BuiltInUnsigned  r[] = { u, v };
Chris@16:             unsigned         which = static_cast<bool>( u & 1u );
Chris@16: 
Chris@16:             // Whittle down the values via their differences
Chris@16:             do
Chris@16:             {
Chris@16: #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
Chris@16:                 while ( !(r[ which ] & 1u) )
Chris@16:                 {
Chris@16:                     r[ which ] = (r[which] >> 1);
Chris@16:                 }
Chris@16: #else
Chris@16:                 // Remove factors of two from the even one
Chris@16:                 while ( !(r[ which ] & 1u) )
Chris@16:                 {
Chris@16:                     r[ which ] >>= 1;
Chris@16:                 }
Chris@16: #endif
Chris@16: 
Chris@16:                 // Replace the larger of the two with their difference
Chris@16:                 if ( r[!which] > r[which] )
Chris@16:                 {
Chris@16:                     which ^= 1u;
Chris@16:                 }
Chris@16: 
Chris@16:                 r[ which ] -= r[ !which ];
Chris@16:             }
Chris@16:             while ( r[which] );
Chris@16: 
Chris@16:             // Shift-in the common factor of 2 to the residues' GCD
Chris@16:             return r[ !which ] << shifts;
Chris@16:         }
Chris@16:         else
Chris@16:         {
Chris@16:             // At least one input is zero, return the other
Chris@16:             // (adding since zero is the additive identity)
Chris@16:             // or zero if both are zero.
Chris@16:             return u + v;
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     // Least common multiple for rings (including unsigned integers)
Chris@16:     template < typename RingType >
Chris@16:     inline
Chris@16:     RingType
Chris@16:     lcm_euclidean
Chris@16:     (
Chris@16:         RingType const &  a,
Chris@16:         RingType const &  b
Chris@16:     )
Chris@16:     {
Chris@16:         RingType const  zero = static_cast<RingType>( 0 );
Chris@16:         RingType const  temp = gcd_euclidean( a, b );
Chris@16: 
Chris@16:         return ( temp != zero ) ? ( a / temp * b ) : zero;
Chris@16:     }
Chris@16: 
Chris@16:     // Least common multiple for (signed) integers
Chris@16:     template < typename IntegerType >
Chris@16:     inline
Chris@16:     IntegerType
Chris@16:     lcm_integer
Chris@16:     (
Chris@16:         IntegerType const &  a,
Chris@16:         IntegerType const &  b
Chris@16:     )
Chris@16:     {
Chris@16:         // Avoid repeated construction
Chris@16:         IntegerType const  zero = static_cast<IntegerType>( 0 );
Chris@16:         IntegerType const  result = lcm_euclidean( a, b );
Chris@16: 
Chris@16:         return ( result < zero ) ? static_cast<IntegerType>(-result) : result;
Chris@16:     }
Chris@16: 
Chris@16:     // Function objects to find the best way of computing GCD or LCM
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:     template < typename T, bool IsSpecialized, bool IsSigned >
Chris@16:     struct gcd_optimal_evaluator_helper_t
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return gcd_euclidean( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     template < typename T >
Chris@16:     struct gcd_optimal_evaluator_helper_t< T, true, true >
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return gcd_integer( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     template < typename T >
Chris@16:     struct gcd_optimal_evaluator
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             typedef ::std::numeric_limits<T>  limits_type;
Chris@16: 
Chris@16:             typedef gcd_optimal_evaluator_helper_t<T,
Chris@16:              limits_type::is_specialized, limits_type::is_signed>  helper_type;
Chris@16: 
Chris@16:             helper_type  solver;
Chris@16: 
Chris@16:             return solver( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: #else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:     template < typename T >
Chris@16:     struct gcd_optimal_evaluator
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return gcd_integer( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: #endif
Chris@16: 
Chris@16:     // Specialize for the built-in integers
Chris@16: #define BOOST_PRIVATE_GCD_UF( Ut )                  \
Chris@16:     template < >  struct gcd_optimal_evaluator<Ut>  \
Chris@16:     {  Ut  operator ()( Ut a, Ut b ) const  { return gcd_binary( a, b ); }  }
Chris@16: 
Chris@16:     BOOST_PRIVATE_GCD_UF( unsigned char );
Chris@16:     BOOST_PRIVATE_GCD_UF( unsigned short );
Chris@16:     BOOST_PRIVATE_GCD_UF( unsigned );
Chris@16:     BOOST_PRIVATE_GCD_UF( unsigned long );
Chris@16: 
Chris@16: #ifdef BOOST_HAS_LONG_LONG
Chris@16:     BOOST_PRIVATE_GCD_UF( boost::ulong_long_type );
Chris@16: #elif defined(BOOST_HAS_MS_INT64)
Chris@16:     BOOST_PRIVATE_GCD_UF( unsigned __int64 );
Chris@16: #endif
Chris@16: 
Chris@16: #if CHAR_MIN == 0
Chris@16:     BOOST_PRIVATE_GCD_UF( char ); // char is unsigned
Chris@16: #endif
Chris@16: 
Chris@16: #undef BOOST_PRIVATE_GCD_UF
Chris@16: 
Chris@16: #define BOOST_PRIVATE_GCD_SF( St, Ut )                            \
Chris@16:     template < >  struct gcd_optimal_evaluator<St>                \
Chris@16:     {  St  operator ()( St a, St b ) const  { Ut const  a_abs =   \
Chris@16:     static_cast<Ut>( a < 0 ? -a : +a ), b_abs = static_cast<Ut>(  \
Chris@16:     b < 0 ? -b : +b ); return static_cast<St>(                    \
Chris@16:     gcd_optimal_evaluator<Ut>()(a_abs, b_abs) ); }  }
Chris@16: 
Chris@16:     BOOST_PRIVATE_GCD_SF( signed char, unsigned char );
Chris@16:     BOOST_PRIVATE_GCD_SF( short, unsigned short );
Chris@16:     BOOST_PRIVATE_GCD_SF( int, unsigned );
Chris@16:     BOOST_PRIVATE_GCD_SF( long, unsigned long );
Chris@16: 
Chris@16: #if CHAR_MIN < 0
Chris@16:     BOOST_PRIVATE_GCD_SF( char, unsigned char ); // char is signed
Chris@16: #endif
Chris@16: 
Chris@16: #ifdef BOOST_HAS_LONG_LONG
Chris@16:     BOOST_PRIVATE_GCD_SF( boost::long_long_type, boost::ulong_long_type );
Chris@16: #elif defined(BOOST_HAS_MS_INT64)
Chris@16:     BOOST_PRIVATE_GCD_SF( __int64, unsigned __int64 );
Chris@16: #endif
Chris@16: 
Chris@16: #undef BOOST_PRIVATE_GCD_SF
Chris@16: 
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:     template < typename T, bool IsSpecialized, bool IsSigned >
Chris@16:     struct lcm_optimal_evaluator_helper_t
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return lcm_euclidean( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     template < typename T >
Chris@16:     struct lcm_optimal_evaluator_helper_t< T, true, true >
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return lcm_integer( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     template < typename T >
Chris@16:     struct lcm_optimal_evaluator
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             typedef ::std::numeric_limits<T>  limits_type;
Chris@16: 
Chris@16:             typedef lcm_optimal_evaluator_helper_t<T,
Chris@16:              limits_type::is_specialized, limits_type::is_signed>  helper_type;
Chris@16: 
Chris@16:             helper_type  solver;
Chris@16: 
Chris@16:             return solver( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: #else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:     template < typename T >
Chris@16:     struct lcm_optimal_evaluator
Chris@16:     {
Chris@16:         T  operator ()( T const &a, T const &b )
Chris@16:         {
Chris@16:             return lcm_integer( a, b );
Chris@16:         }
Chris@16:     };
Chris@16: #endif
Chris@16: 
Chris@16:     // Functions to find the GCD or LCM in the best way
Chris@16:     template < typename T >
Chris@16:     inline
Chris@16:     T
Chris@16:     gcd_optimal
Chris@16:     (
Chris@16:         T const &  a,
Chris@16:         T const &  b
Chris@16:     )
Chris@16:     {
Chris@16:         gcd_optimal_evaluator<T>  solver;
Chris@16: 
Chris@16:         return solver( a, b );
Chris@16:     }
Chris@16: 
Chris@16:     template < typename T >
Chris@16:     inline
Chris@16:     T
Chris@16:     lcm_optimal
Chris@16:     (
Chris@16:         T const &  a,
Chris@16:         T const &  b
Chris@16:     )
Chris@16:     {
Chris@16:         lcm_optimal_evaluator<T>  solver;
Chris@16: 
Chris@16:         return solver( a, b );
Chris@16:     }
Chris@16: 
Chris@16: }  // namespace detail
Chris@16: 
Chris@16: 
Chris@16: //  Greatest common divisor evaluator member function definition  ------------//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: inline
Chris@16: typename gcd_evaluator<IntegerType>::result_type
Chris@16: gcd_evaluator<IntegerType>::operator ()
Chris@16: (
Chris@16:     first_argument_type const &   a,
Chris@16:     second_argument_type const &  b
Chris@16: ) const
Chris@16: {
Chris@16:     return detail::gcd_optimal( a, b );
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: //  Least common multiple evaluator member function definition  --------------//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: inline
Chris@16: typename lcm_evaluator<IntegerType>::result_type
Chris@16: lcm_evaluator<IntegerType>::operator ()
Chris@16: (
Chris@16:     first_argument_type const &   a,
Chris@16:     second_argument_type const &  b
Chris@16: ) const
Chris@16: {
Chris@16:     return detail::lcm_optimal( a, b );
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: //  Greatest common divisor and least common multiple function definitions  --//
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: inline
Chris@16: IntegerType
Chris@16: gcd
Chris@16: (
Chris@16:     IntegerType const &  a,
Chris@16:     IntegerType const &  b
Chris@16: )
Chris@16: {
Chris@16:     gcd_evaluator<IntegerType>  solver;
Chris@16: 
Chris@16:     return solver( a, b );
Chris@16: }
Chris@16: 
Chris@16: template < typename IntegerType >
Chris@16: inline
Chris@16: IntegerType
Chris@16: lcm
Chris@16: (
Chris@16:     IntegerType const &  a,
Chris@16:     IntegerType const &  b
Chris@16: )
Chris@16: {
Chris@16:     lcm_evaluator<IntegerType>  solver;
Chris@16: 
Chris@16:     return solver( a, b );
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: }  // namespace math
Chris@16: }  // namespace boost
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: 
Chris@16: #endif  // BOOST_MATH_COMMON_FACTOR_RT_HPP