boost/integer/mod_inverse.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_MOD_INVERSE_HPP
0008 #define BOOST_INTEGER_MOD_INVERSE_HPP
0009 #include <stdexcept>
0010 #include <boost/throw_exception.hpp>
0011 #include <boost/integer/extended_euclidean.hpp>
0012
0013 namespace boost { namespace integer {
0014
0015 // From "The Joy of Factoring", Algorithm 2.7.
0016 // Here's some others names I've found for this function:
0017 // PowerMod[a, -1, m] (Mathematica)
0018 // mpz_invert (gmplib)
0019 // modinv (some dude on stackoverflow)
0020 // Would mod_inverse be sometimes mistaken as the modular *additive* inverse?
0021 // In any case, I think this is the best name we can get for this function without agonizing.
0022 template<class Z>
0023 Z mod_inverse(Z a, Z modulus)
0024 {
0025     if (modulus < Z(2))
0026     {
0027         BOOST_THROW_EXCEPTION(std::domain_error("mod_inverse: modulus must be > 1"));
0028     }
0029     // make sure a < modulus:
0030     a = a % modulus;
0031     if (a == Z(0))
0032     {
0033         // a doesn't have a modular multiplicative inverse:
0034         return Z(0);
0035     }
0036     boost::integer::euclidean_result_t<Z> u = boost::integer::extended_euclidean(a, modulus);
0037     if (u.gcd > Z(1))
0038     {
0039         return Z(0);
0040     }
0041     // x might not be in the range 0 < x < m, let's fix that:
0042     while (u.x <= Z(0))
0043     {
0044         u.x += modulus;
0045     }
0046     // While indeed this is an inexpensive and comforting check,
0047     // the multiplication overflows and hence makes the check itself buggy.
0048     //BOOST_ASSERT(u.x*a % modulus == 1);
0049     return u.x;
0050 }
0051
0052 }}
0053 #endif