common_factor_rt.hpp 13 KB

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