mod_inverse.hpp 1.6 KB

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