boost/math/common_factor_rt.hpp

   1 //  Boost common_factor_rt.hpp header file  ----------------------------------//
   2
   3 //  (C) Copyright Daryle Walker and Paul Moore 2001-2002.  Permission to copy,
   4 //  use, modify, sell and distribute this software is granted provided this
   5 //  copyright notice appears in all copies.  This software is provided "as is"
   6 //  without express or implied warranty, and with no claim as to its suitability
   7 //  for any purpose.
   8
   9 // boostinspect:nolicense (don't complain about the lack of a Boost license)
  10 // (Paul Moore hasn't been in contact for years, so there's no way to change the
  11 // license.)
  12
  13 //  See http://www.boost.org for updates, documentation, and revision history.
  14
  15 #ifndef BOOST_MATH_COMMON_FACTOR_RT_HPP
  16 #define BOOST_MATH_COMMON_FACTOR_RT_HPP
  17
  18 #include <boost/math_fwd.hpp>  // self include
  19
  20 #include <boost/config.hpp>  // for BOOST_NESTED_TEMPLATE, etc.
  21 #include <boost/limits.hpp>  // for std::numeric_limits
  22 #include <climits>           // for CHAR_MIN
  23 #include <boost/detail/workaround.hpp>
  24
  25 #ifdef BOOST_MSVC
  26 #pragma warning(push)
  27 #pragma warning(disable:4127 4244)  // Conditional expression is constant
  28 #endif
  29
  30 namespace boost
  31 {
  32 namespace math
  33 {
  34
  35
  36 //  Forward declarations for function templates  -----------------------------//
  37
  38 template < typename IntegerType >
  39     IntegerType  gcd( IntegerType const &a, IntegerType const &b );
  40
  41 template < typename IntegerType >
  42     IntegerType  lcm( IntegerType const &a, IntegerType const &b );
  43
  44
  45 //  Greatest common divisor evaluator class declaration  ---------------------//
  46
  47 template < typename IntegerType >
  48 class gcd_evaluator
  49 {
  50 public:
  51     // Types
  52     typedef IntegerType  result_type, first_argument_type, second_argument_type;
  53
  54     // Function object interface
  55     result_type  operator ()( first_argument_type const &a,
  56      second_argument_type const &b ) const;
  57
  58 };  // boost::math::gcd_evaluator
  59
  60
  61 //  Least common multiple evaluator class declaration  -----------------------//
  62
  63 template < typename IntegerType >
  64 class lcm_evaluator
  65 {
  66 public:
  67     // Types
  68     typedef IntegerType  result_type, first_argument_type, second_argument_type;
  69
  70     // Function object interface
  71     result_type  operator ()( first_argument_type const &a,
  72      second_argument_type const &b ) const;
  73
  74 };  // boost::math::lcm_evaluator
  75
  76
  77 //  Implementation details  --------------------------------------------------//
  78
  79 namespace detail
  80 {
  81     // Greatest common divisor for rings (including unsigned integers)
  82     template < typename RingType >
  83     RingType
  84     gcd_euclidean
  85     (
  86         RingType a,
  87         RingType b
  88     )
  89     {
  90         // Avoid repeated construction
  91         #ifndef __BORLANDC__
  92         RingType const  zero = static_cast<RingType>( 0 );
  93         #else
  94         RingType  zero = static_cast<RingType>( 0 );
  95         #endif
  96
  97         // Reduce by GCD-remainder property [GCD(a,b) == GCD(b,a MOD b)]
  98         while ( true )
  99         {
 100             if ( a == zero )
 101                 return b;
 102             b %= a;
 103
 104             if ( b == zero )
 105                 return a;
 106             a %= b;
 107         }
 108     }
 109
 110     // Greatest common divisor for (signed) integers
 111     template < typename IntegerType >
 112     inline
 113     IntegerType
 114     gcd_integer
 115     (
 116         IntegerType const &  a,
 117         IntegerType const &  b
 118     )
 119     {
 120         // Avoid repeated construction
 121         IntegerType const  zero = static_cast<IntegerType>( 0 );
 122         IntegerType const  result = gcd_euclidean( a, b );
 123
 124         return ( result < zero ) ? static_cast<IntegerType>(-result) : result;
 125     }
 126
 127     // Greatest common divisor for unsigned binary integers
 128     template < typename BuiltInUnsigned >
 129     BuiltInUnsigned
 130     gcd_binary
 131     (
 132         BuiltInUnsigned  u,
 133         BuiltInUnsigned  v
 134     )
 135     {
 136         if ( u && v )
 137         {
 138             // Shift out common factors of 2
 139             unsigned  shifts = 0;
 140
 141             while ( !(u & 1u) && !(v & 1u) )
 142             {
 143                 ++shifts;
 144                 u >>= 1;
 145                 v >>= 1;
 146             }
 147
 148             // Start with the still-even one, if any
 149             BuiltInUnsigned  r[] = { u, v };
 150             unsigned         which = static_cast<bool>( u & 1u );
 151
 152             // Whittle down the values via their differences
 153             do
 154             {
 155 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 156                 while ( !(r[ which ] & 1u) )
 157                 {
 158                     r[ which ] = (r[which] >> 1);
 159                 }
 160 #else
 161                 // Remove factors of two from the even one
 162                 while ( !(r[ which ] & 1u) )
 163                 {
 164                     r[ which ] >>= 1;
 165                 }
 166 #endif
 167
 168                 // Replace the larger of the two with their difference
 169                 if ( r[!which] > r[which] )
 170                 {
 171                     which ^= 1u;
 172                 }
 173
 174                 r[ which ] -= r[ !which ];
 175             }
 176             while ( r[which] );
 177
 178             // Shift-in the common factor of 2 to the residues' GCD
 179             return r[ !which ] << shifts;
 180         }
 181         else
 182         {
 183             // At least one input is zero, return the other
 184             // (adding since zero is the additive identity)
 185             // or zero if both are zero.
 186             return u + v;
 187         }
 188     }
 189
 190     // Least common multiple for rings (including unsigned integers)
 191     template < typename RingType >
 192     inline
 193     RingType
 194     lcm_euclidean
 195     (
 196         RingType const &  a,
 197         RingType const &  b
 198     )
 199     {
 200         RingType const  zero = static_cast<RingType>( 0 );
 201         RingType const  temp = gcd_euclidean( a, b );
 202
 203         return ( temp != zero ) ? ( a / temp * b ) : zero;
 204     }
 205
 206     // Least common multiple for (signed) integers
 207     template < typename IntegerType >
 208     inline
 209     IntegerType
 210     lcm_integer
 211     (
 212         IntegerType const &  a,
 213         IntegerType const &  b
 214     )
 215     {
 216         // Avoid repeated construction
 217         IntegerType const  zero = static_cast<IntegerType>( 0 );
 218         IntegerType const  result = lcm_euclidean( a, b );
 219
 220         return ( result < zero ) ? static_cast<IntegerType>(-result) : result;
 221     }
 222
 223     // Function objects to find the best way of computing GCD or LCM
 224 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 225     template < typename T, bool IsSpecialized, bool IsSigned >
 226     struct gcd_optimal_evaluator_helper_t
 227     {
 228         T  operator ()( T const &a, T const &b )
 229         {
 230             return gcd_euclidean( a, b );
 231         }
 232     };
 233
 234     template < typename T >
 235     struct gcd_optimal_evaluator_helper_t< T, true, true >
 236     {
 237         T  operator ()( T const &a, T const &b )
 238         {
 239             return gcd_integer( a, b );
 240         }
 241     };
 242
 243     template < typename T >
 244     struct gcd_optimal_evaluator
 245     {
 246         T  operator ()( T const &a, T const &b )
 247         {
 248             typedef ::std::numeric_limits<T>  limits_type;
 249
 250             typedef gcd_optimal_evaluator_helper_t<T,
 251              limits_type::is_specialized, limits_type::is_signed>  helper_type;
 252
 253             helper_type  solver;
 254
 255             return solver( a, b );
 256         }
 257     };
 258 #else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 259     template < typename T >
 260     struct gcd_optimal_evaluator
 261     {
 262         T  operator ()( T const &a, T const &b )
 263         {
 264             return gcd_integer( a, b );
 265         }
 266     };
 267 #endif
 268
 269     // Specialize for the built-in integers
 270 #define BOOST_PRIVATE_GCD_UF( Ut )                  \
 271     template < >  struct gcd_optimal_evaluator<Ut>  \
 272     {  Ut  operator ()( Ut a, Ut b ) const  { return gcd_binary( a, b ); }  }
 273
 274     BOOST_PRIVATE_GCD_UF( unsigned char );
 275     BOOST_PRIVATE_GCD_UF( unsigned short );
 276     BOOST_PRIVATE_GCD_UF( unsigned );
 277     BOOST_PRIVATE_GCD_UF( unsigned long );
 278
 279 #ifdef BOOST_HAS_LONG_LONG
 280     BOOST_PRIVATE_GCD_UF( boost::ulong_long_type );
 281 #elif defined(BOOST_HAS_MS_INT64)
 282     BOOST_PRIVATE_GCD_UF( unsigned __int64 );
 283 #endif
 284
 285 #if CHAR_MIN == 0
 286     BOOST_PRIVATE_GCD_UF( char ); // char is unsigned
 287 #endif
 288
 289 #undef BOOST_PRIVATE_GCD_UF
 290
 291 #define BOOST_PRIVATE_GCD_SF( St, Ut )                            \
 292     template < >  struct gcd_optimal_evaluator<St>                \
 293     {  St  operator ()( St a, St b ) const  { Ut const  a_abs =   \
 294     static_cast<Ut>( a < 0 ? -a : +a ), b_abs = static_cast<Ut>(  \
 295     b < 0 ? -b : +b ); return static_cast<St>(                    \
 296     gcd_optimal_evaluator<Ut>()(a_abs, b_abs) ); }  }
 297
 298     BOOST_PRIVATE_GCD_SF( signed char, unsigned char );
 299     BOOST_PRIVATE_GCD_SF( short, unsigned short );
 300     BOOST_PRIVATE_GCD_SF( int, unsigned );
 301     BOOST_PRIVATE_GCD_SF( long, unsigned long );
 302
 303 #if CHAR_MIN < 0
 304     BOOST_PRIVATE_GCD_SF( char, unsigned char ); // char is signed
 305 #endif
 306
 307 #ifdef BOOST_HAS_LONG_LONG
 308     BOOST_PRIVATE_GCD_SF( boost::long_long_type, boost::ulong_long_type );
 309 #elif defined(BOOST_HAS_MS_INT64)
 310     BOOST_PRIVATE_GCD_SF( __int64, unsigned __int64 );
 311 #endif
 312
 313 #undef BOOST_PRIVATE_GCD_SF
 314
 315 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 316     template < typename T, bool IsSpecialized, bool IsSigned >
 317     struct lcm_optimal_evaluator_helper_t
 318     {
 319         T  operator ()( T const &a, T const &b )
 320         {
 321             return lcm_euclidean( a, b );
 322         }
 323     };
 324
 325     template < typename T >
 326     struct lcm_optimal_evaluator_helper_t< T, true, true >
 327     {
 328         T  operator ()( T const &a, T const &b )
 329         {
 330             return lcm_integer( a, b );
 331         }
 332     };
 333
 334     template < typename T >
 335     struct lcm_optimal_evaluator
 336     {
 337         T  operator ()( T const &a, T const &b )
 338         {
 339             typedef ::std::numeric_limits<T>  limits_type;
 340
 341             typedef lcm_optimal_evaluator_helper_t<T,
 342              limits_type::is_specialized, limits_type::is_signed>  helper_type;
 343
 344             helper_type  solver;
 345
 346             return solver( a, b );
 347         }
 348     };
 349 #else // BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 350     template < typename T >
 351     struct lcm_optimal_evaluator
 352     {
 353         T  operator ()( T const &a, T const &b )
 354         {
 355             return lcm_integer( a, b );
 356         }
 357     };
 358 #endif
 359
 360     // Functions to find the GCD or LCM in the best way
 361     template < typename T >
 362     inline
 363     T
 364     gcd_optimal
 365     (
 366         T const &  a,
 367         T const &  b
 368     )
 369     {
 370         gcd_optimal_evaluator<T>  solver;
 371
 372         return solver( a, b );
 373     }
 374
 375     template < typename T >
 376     inline
 377     T
 378     lcm_optimal
 379     (
 380         T const &  a,
 381         T const &  b
 382     )
 383     {
 384         lcm_optimal_evaluator<T>  solver;
 385
 386         return solver( a, b );
 387     }
 388
 389 }  // namespace detail
 390
 391
 392 //  Greatest common divisor evaluator member function definition  ------------//
 393
 394 template < typename IntegerType >
 395 inline
 396 typename gcd_evaluator<IntegerType>::result_type
 397 gcd_evaluator<IntegerType>::operator ()
 398 (
 399     first_argument_type const &   a,
 400     second_argument_type const &  b
 401 ) const
 402 {
 403     return detail::gcd_optimal( a, b );
 404 }
 405
 406
 407 //  Least common multiple evaluator member function definition  --------------//
 408
 409 template < typename IntegerType >
 410 inline
 411 typename lcm_evaluator<IntegerType>::result_type
 412 lcm_evaluator<IntegerType>::operator ()
 413 (
 414     first_argument_type const &   a,
 415     second_argument_type const &  b
 416 ) const
 417 {
 418     return detail::lcm_optimal( a, b );
 419 }
 420
 421
 422 //  Greatest common divisor and least common multiple function definitions  --//
 423
 424 template < typename IntegerType >
 425 inline
 426 IntegerType
 427 gcd
 428 (
 429     IntegerType const &  a,
 430     IntegerType const &  b
 431 )
 432 {
 433     gcd_evaluator<IntegerType>  solver;
 434
 435     return solver( a, b );
 436 }
 437
 438 template < typename IntegerType >
 439 inline
 440 IntegerType
 441 lcm
 442 (
 443     IntegerType const &  a,
 444     IntegerType const &  b
 445 )
 446 {
 447     lcm_evaluator<IntegerType>  solver;
 448
 449     return solver( a, b );
 450 }
 451
 452
 453 }  // namespace math
 454 }  // namespace boost
 455
 456 #ifdef BOOST_MSVC
 457 #pragma warning(pop)
 458 #endif
 459
 460 #endif  // BOOST_MATH_COMMON_FACTOR_RT_HPP