diff --git a/include/boost/integer/extended_euclidean.hpp b/include/boost/integer/extended_euclidean.hpp index 64d10d9..93c8c54 100644 --- a/include/boost/integer/extended_euclidean.hpp +++ b/include/boost/integer/extended_euclidean.hpp @@ -13,9 +13,8 @@ namespace boost { namespace integer { -// From "The Joy of Factoring", Algorithm 2.7. +// From "The Joy of Factoring", Algorithm 2.7, with a small optimization to remove tmps from Wikipedia. // Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y). -// Is this the natural ordering?, or must people simply have to read the docs? template struct euclidean_result_t { @@ -38,45 +37,27 @@ euclidean_result_t extended_euclidean(Z m, Z n) { BOOST_THROW_EXCEPTION(std::domain_error("Arguments must be strictly positive.\n")); } - bool swapped = false; - if (m < n) - { - swapped = true; - std::swap(m, n); - } - Z u0 = m; - Z u1 = 1; - Z u2 = 0; - Z v0 = n; - Z v1 = 0; - Z v2 = 1; - Z w0; - Z w1; - Z w2; - while(v0 > 0) - { - Z q = u0/v0; - w0 = u0 - q*v0; - w1 = u1 - q*v1; - w2 = u2 - q*v2; - u0 = v0; - u1 = v1; - u2 = v2; - v0 = w0; - v1 = w1; - v2 = w2; - } - if (swapped) - { - std::swap(u1, u2); - BOOST_ASSERT(u2*m+u1*n==u0); - } - else - { - BOOST_ASSERT(u1*m+u2*n==u0); + + Z s = 0; + Z old_s = 1; + Z r = n; + Z old_r = m; + + while (r != 0) { + Z q = old_r/r; + Z tmp = r; + r = old_r - q*tmp; + old_r = tmp; + + tmp = s; + s = old_s - q*tmp; + old_s = tmp; } - return {u0, u1, u2}; + Z y = (old_r - old_s*m)/n; + + BOOST_ASSERT(old_s*m+y*n==old_r); + return {old_r, old_s, y}; } }}