boost/integer/extended_euclidean.hpp

0001 /*
0002  *  (C) Copyright Nick Thompson 2018.
0003  *  Use, modification and distribution are subject to the
0004  *  Boost Software License, Version 1.0. (See accompanying file
0005  *  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
0006  */
0007 #ifndef BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
0008 #define BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
0009 #include <limits>
0010 #include <stdexcept>
0011 #include <boost/throw_exception.hpp>
0012 #include <boost/core/invoke_swap.hpp>
0013 #include <boost/core/enable_if.hpp>
0014
0015 namespace boost { namespace integer {
0016
0017 // From "The Joy of Factoring", Algorithm 2.7, with a small optimization to remove tmps from Wikipedia.
0018 // Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y).
0019
0020 template<class Z>
0021 struct euclidean_result_t
0022 {
0023     Z gcd;
0024     Z x;
0025     Z y;
0026 };
0027
0028 template<class Z>
0029 typename boost::enable_if_c< std::numeric_limits< Z >::is_signed, euclidean_result_t< Z > >::type
0030 extended_euclidean(Z m, Z n)
0031 {
0032     if (m < 1 || n < 1)
0033     {
0034         BOOST_THROW_EXCEPTION(std::domain_error("extended_euclidean: arguments must be strictly positive"));
0035     }
0036
0037     bool swapped = false;
0038     if (m < n)
0039     {
0040         swapped = true;
0041         boost::core::invoke_swap(m, n);
0042     }
0043     Z u0 = m;
0044     Z u1 = 1;
0045     Z u2 = 0;
0046     Z v0 = n;
0047     Z v1 = 0;
0048     Z v2 = 1;
0049     Z w0;
0050     Z w1;
0051     Z w2;
0052     while(v0 > 0)
0053     {
0054         Z q = u0/v0;
0055         w0 = u0 - q*v0;
0056         w1 = u1 - q*v1;
0057         w2 = u2 - q*v2;
0058         u0 = v0;
0059         u1 = v1;
0060         u2 = v2;
0061         v0 = w0;
0062         v1 = w1;
0063         v2 = w2;
0064     }
0065
0066     euclidean_result_t< Z > result;
0067     result.gcd = u0;
0068     if (!swapped)
0069     {
0070         result.x = u1;
0071         result.y = u2;
0072     }
0073     else
0074     {
0075         result.x = u2;
0076         result.y = u1;
0077     }
0078
0079     return result;
0080 }
0081
0082 }}
0083 #endif