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_MOD_INVERSE_HPP
8 #define BOOST_INTEGER_MOD_INVERSE_HPP
9 #include <stdexcept>
10 #include <boost/throw_exception.hpp>
11 #include <boost/integer/extended_euclidean.hpp>
12
13 namespace boost { namespace integer {
14
15 // From "The Joy of Factoring", Algorithm 2.7.
16 // Here's some others names I've found for this function:
17 // PowerMod[a, -1, m] (Mathematica)
18 // mpz_invert (gmplib)
19 // modinv (some dude on stackoverflow)
20 // Would mod_inverse be sometimes mistaken as the modular *additive* inverse?
21 // In any case, I think this is the best name we can get for this function without agonizing.
22 template<class Z>
mod_inverse(Z a,Z modulus)23 Z mod_inverse(Z a, Z modulus)
24 {
25 if (modulus < Z(2))
26 {
27 BOOST_THROW_EXCEPTION(std::domain_error("mod_inverse: modulus must be > 1"));
28 }
29 // make sure a < modulus:
30 a = a % modulus;
31 if (a == Z(0))
32 {
33 // a doesn't have a modular multiplicative inverse:
34 return Z(0);
35 }
36 boost::integer::euclidean_result_t<Z> u = boost::integer::extended_euclidean(a, modulus);
37 if (u.gcd > Z(1))
38 {
39 return Z(0);
40 }
41 // x might not be in the range 0 < x < m, let's fix that:
42 while (u.x <= Z(0))
43 {
44 u.x += modulus;
45 }
46 // While indeed this is an inexpensive and comforting check,
47 // the multiplication overflows and hence makes the check itself buggy.
48 //BOOST_ASSERT(u.x*a % modulus == 1);
49 return u.x;
50 }
51
52 }}
53 #endif
54