extended_euclidean.hpp 1.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283
  1. /*
  2. * (C) Copyright Nick Thompson 2018.
  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 https://www.boost.org/LICENSE_1_0.txt)
  6. */
  7. #ifndef BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
  8. #define BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
  9. #include <limits>
  10. #include <stdexcept>
  11. #include <boost/throw_exception.hpp>
  12. #include <boost/core/swap.hpp>
  13. #include <boost/core/enable_if.hpp>
  14. namespace boost { namespace integer {
  15. // From "The Joy of Factoring", Algorithm 2.7, with a small optimization to remove tmps from Wikipedia.
  16. // Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y).
  17. template<class Z>
  18. struct euclidean_result_t
  19. {
  20. Z gcd;
  21. Z x;
  22. Z y;
  23. };
  24. template<class Z>
  25. typename boost::enable_if_c< std::numeric_limits< Z >::is_signed, euclidean_result_t< Z > >::type
  26. extended_euclidean(Z m, Z n)
  27. {
  28. if (m < 1 || n < 1)
  29. {
  30. BOOST_THROW_EXCEPTION(std::domain_error("extended_euclidean: arguments must be strictly positive"));
  31. }
  32. bool swapped = false;
  33. if (m < n)
  34. {
  35. swapped = true;
  36. boost::swap(m, n);
  37. }
  38. Z u0 = m;
  39. Z u1 = 1;
  40. Z u2 = 0;
  41. Z v0 = n;
  42. Z v1 = 0;
  43. Z v2 = 1;
  44. Z w0;
  45. Z w1;
  46. Z w2;
  47. while(v0 > 0)
  48. {
  49. Z q = u0/v0;
  50. w0 = u0 - q*v0;
  51. w1 = u1 - q*v1;
  52. w2 = u2 - q*v2;
  53. u0 = v0;
  54. u1 = v1;
  55. u2 = v2;
  56. v0 = w0;
  57. v1 = w1;
  58. v2 = w2;
  59. }
  60. euclidean_result_t< Z > result;
  61. result.gcd = u0;
  62. if (!swapped)
  63. {
  64. result.x = u1;
  65. result.y = u2;
  66. }
  67. else
  68. {
  69. result.x = u2;
  70. result.y = u1;
  71. }
  72. return result;
  73. }
  74. }}
  75. #endif