]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 http://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 | ||
15 | namespace boost { namespace integer { | |
16 | ||
17 | // From "The Joy of Factoring", Algorithm 2.7, with a small optimization to remove tmps from Wikipedia. | |
18 | // Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y). | |
19 | ||
20 | template<class Z> | |
21 | struct euclidean_result_t | |
22 | { | |
23 | Z gcd; | |
24 | Z x; | |
25 | Z y; | |
26 | }; | |
27 | ||
28 | template<class Z> | |
29 | typename boost::enable_if_c< std::numeric_limits< Z >::is_signed, euclidean_result_t< Z > >::type | |
30 | extended_euclidean(Z m, Z n) | |
31 | { | |
32 | if (m < 1 || n < 1) | |
33 | { | |
34 | BOOST_THROW_EXCEPTION(std::domain_error("extended_euclidean: arguments must be strictly positive")); | |
35 | } | |
36 | ||
37 | bool swapped = false; | |
38 | if (m < n) | |
39 | { | |
40 | swapped = true; | |
41 | boost::swap(m, n); | |
42 | } | |
43 | Z u0 = m; | |
44 | Z u1 = 1; | |
45 | Z u2 = 0; | |
46 | Z v0 = n; | |
47 | Z v1 = 0; | |
48 | Z v2 = 1; | |
49 | Z w0; | |
50 | Z w1; | |
51 | Z w2; | |
52 | while(v0 > 0) | |
53 | { | |
54 | Z q = u0/v0; | |
55 | w0 = u0 - q*v0; | |
56 | w1 = u1 - q*v1; | |
57 | w2 = u2 - q*v2; | |
58 | u0 = v0; | |
59 | u1 = v1; | |
60 | u2 = v2; | |
61 | v0 = w0; | |
62 | v1 = w1; | |
63 | v2 = w2; | |
64 | } | |
65 | ||
66 | euclidean_result_t< Z > result; | |
67 | result.gcd = u0; | |
68 | if (!swapped) | |
69 | { | |
70 | result.x = u1; | |
71 | result.y = u2; | |
72 | } | |
73 | else | |
74 | { | |
75 | result.x = u2; | |
76 | result.y = u1; | |
77 | } | |
78 | ||
79 | return result; | |
80 | } | |
81 | ||
82 | }} | |
83 | #endif |