Modular multiplicative inverse function in Python
Asked Answered
H

13

189

Does some standard Python module contain a function to compute modular multiplicative inverse of a number, i.e. a number y = invmod(x, p) such that x*y == 1 (mod p)? Google doesn't seem to give any good hints on this.

Of course, one can come up with home-brewed 10-liner of extended Euclidean algorithm, but why reinvent the wheel.

For example, Java's BigInteger has modInverse method. Doesn't Python have something similar?

Hammered answered 25/1, 2011 at 20:50 Comment(2)
In Python 3.8 (due to be released later this year), you'll be able to use the built-in pow function for this: y = pow(x, -1, p). See bugs.python.org/issue36027. It only took 8.5 years from the question being asked to a solution appearing in the standard library!Velvety
I see @MarkDickinson modestly neglected to mention that ey is the author of this very useful enhancement, so I will. Thanks for this work, Mark, it looks great!Variate
H
240

Python 3.8+

y = pow(x, -1, p)

Python 3.7 and earlier

Maybe someone will find this useful (from wikibooks):

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m
Homologue answered 18/3, 2012 at 12:8 Comment(4)
I was having problems with negative numbers using this algorithm. modinv(-3, 11) didn't work. I fixed it by replacing egcd with the implementation on page two of this pdf: anh.cs.luc.edu/331/notes/xgcd.pdf Hope that helps!Johnston
@Johnston You can also just reduce -3 modulo 11 to make it positive, in this case modinv(-3, 11) == modinv(-3 + 11, 11) == modinv(8, 11). That's probably what the algorithm in your PDF happens to do at some point.Glanti
If you happen to be using sympy, then x, _, g = sympy.numbers.igcdex(a, m) does the trick.Rattrap
Love that it's baked into the language for 3.8+Leister
D
67

If your modulus is prime (you call it p) then you may simply compute:

y = x**(p-2) mod p  # Pseudocode

Or in Python proper:

y = pow(x, p-2, p)

Here is someone who has implemented some number theory capabilities in Python: http://www.math.umbc.edu/~campbell/Computers/Python/numbthy.html

Here is an example done at the prompt:

m = 1000000007
x = 1234567
y = pow(x,m-2,m)
y
989145189L
x*y
1221166008548163L
x*y % m
1L
Digressive answered 25/1, 2011 at 21:2 Comment(5)
Naive exponentiation is not an option because of time (and memory) limit for any reasonably big value of p like say 1000000007.Hammered
modular exponentiation is done with at most N*2 multiplications where N is the number of bits in the exponent. using a modulus of 2**63-1 the inverse can be computed at the prompt and returns a result immediately.Digressive
Wow, awesome. I'm aware of quick exponentiation, I just wasn't aware that pow() function can take third argument which turns it into modular exponentiation.Hammered
for large primes, using the egcd method in python is actually faster than using pow(a, m-2, m) ....Beehive
By the way this works because from Fermat little theorem pow(x,m-1,m) must be 1. Hence (pow(x,m-2,m) * x) % m == 1. So pow(x,m-2,m) is the inverse of x (mod m).Preengage
B
25

You might also want to look at the gmpy module. It is an interface between Python and the GMP multiple-precision library. gmpy provides an invert function that does exactly what you need:

>>> import gmpy
>>> gmpy.invert(1234567, 1000000007)
mpz(989145189)

Updated answer

As noted by @hyh , the gmpy.invert() returns 0 if the inverse does not exist. That matches the behavior of GMP's mpz_invert() function. gmpy.divm(a, b, m) provides a general solution to a=bx (mod m).

>>> gmpy.divm(1, 1234567, 1000000007)
mpz(989145189)
>>> gmpy.divm(1, 0, 5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: not invertible
>>> gmpy.divm(1, 4, 8)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: not invertible
>>> gmpy.divm(1, 4, 9)
mpz(7)

divm() will return a solution when gcd(b,m) == 1 and raises an exception when the multiplicative inverse does not exist.

Disclaimer: I'm the current maintainer of the gmpy library.

Updated answer 2

gmpy2 now properly raises an exception when the inverse does not exists:

>>> import gmpy2

>>> gmpy2.invert(0,5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: invert() no inverse exists
Blind answered 26/1, 2011 at 4:13 Comment(7)
This is cool until I found gmpy.invert(0,5) = mpz(0) instead of raising an error...Camm
@hyh Can you report this as an issue at gmpy's home page? It's always appreciated if issues are reported.Blind
BTW, is there a modular multiplication in this gmpy package? (i.e. some function that has the same value but is faster than (a * b)% p ?)Camm
It has been proposed before and I am experimenting with different methods. The simplest approach of just calculating (a * b) % p in a function isn't faster than just evaluating (a * b) % p in Python. The overhead for a function call is greater than the cost of evaluating the expression. See code.google.com/p/gmpy/issues/detail?id=61 for more details.Blind
I see. (Actually in lots of situation, e.g. finite field arithmetic, the modulos is fixed. I hope there is a way to set an environment variable, p such that all arithmetic is then done modulo p. Maybe that needs some bottom up design instead of a few lines of scripts.)Camm
That is one of the options. Actually, I'm working on a re-design to provide better support for threading. After that is done, it will be easy to add a modulus option to context.Blind
The great thing is that this also works for non-prime modulii.Quadratics
C
21

As of 3.8 pythons pow() function can take a modulus and a negative integer. See here. Their case for how to use it is

>>> pow(38, -1, 97)
23
>>> 23 * 38 % 97 == 1
True
Cloudcapped answered 6/11, 2019 at 14:44 Comment(0)
H
10

Here is a one-liner for CodeFights; it is one of the shortest solutions:

MMI = lambda A, n,s=1,t=0,N=0: (n < 2 and t%N or MMI(n, A%n, t, s-A//n*t, N or n),-1)[n<1]

It will return -1 if A has no multiplicative inverse in n.

Usage:

MMI(23, 99) # returns 56
MMI(18, 24) # return -1

The solution uses the Extended Euclidean Algorithm.

Hostess answered 21/4, 2015 at 3:3 Comment(0)
C
9

Sympy, a python module for symbolic mathematics, has a built-in modular inverse function if you don't want to implement your own (or if you're using Sympy already):

from sympy import mod_inverse

mod_inverse(11, 35) # returns 16
mod_inverse(15, 35) # raises ValueError: 'inverse of 15 (mod 35) does not exist'

This doesn't seem to be documented on the Sympy website, but here's the docstring: Sympy mod_inverse docstring on Github

Clamorous answered 8/4, 2018 at 12:44 Comment(1)
mod_inverse is documented on sympy's website: docs.sympy.org/latest/modules/…Charcoal
V
5

Here is a concise 1-liner that does it, without using any external libraries.

# Given 0<a<b, returns the unique c such that 0<c<b and a*c == gcd(a,b) (mod b).
# In particular, if a,b are relatively prime, returns the inverse of a modulo b.
def invmod(a,b): return 0 if a==0 else 1 if b%a==0 else b - invmod(b%a,a)*b//a

Note that this is really just egcd, streamlined to return only the single coefficient of interest.

Variate answered 24/7, 2019 at 1:10 Comment(0)
D
3

I try different solutions from this thread and in the end I use this one:

def egcd(a, b):
    lastremainder, remainder = abs(a), abs(b)
    x, lastx, y, lasty = 0, 1, 1, 0
    while remainder:
        lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
        x, lastx = lastx - quotient*x, x
        y, lasty = lasty - quotient*y, y
    return lastremainder, lastx * (-1 if a < 0 else 1), lasty * (-1 if b < 0 else 1)


def modinv(a, m):
    g, x, y = self.egcd(a, m)
    if g != 1:
        raise ValueError('modinv for {} does not exist'.format(a))
    return x % m

Modular_inverse in Python

Dorkus answered 5/2, 2019 at 12:24 Comment(0)
K
2

Here is my code, it might be sloppy but it seems to work for me anyway.

# a is the number you want the inverse for
# b is the modulus

def mod_inverse(a, b):
    r = -1
    B = b
    A = a
    eq_set = []
    full_set = []
    mod_set = []

    #euclid's algorithm
    while r!=1 and r!=0:
        r = b%a
        q = b//a
        eq_set = [r, b, a, q*-1]
        b = a
        a = r
        full_set.append(eq_set)

    for i in range(0, 4):
        mod_set.append(full_set[-1][i])

    mod_set.insert(2, 1)
    counter = 0

    #extended euclid's algorithm
    for i in range(1, len(full_set)):
        if counter%2 == 0:
            mod_set[2] = full_set[-1*(i+1)][3]*mod_set[4]+mod_set[2]
            mod_set[3] = full_set[-1*(i+1)][1]

        elif counter%2 != 0:
            mod_set[4] = full_set[-1*(i+1)][3]*mod_set[2]+mod_set[4]
            mod_set[1] = full_set[-1*(i+1)][1]

        counter += 1

    if mod_set[3] == B:
        return mod_set[2]%B
    return mod_set[4]%B
Knuckle answered 21/2, 2015 at 0:58 Comment(0)
H
2

The code above will not run in python3 and is less efficient compared to the GCD variants. However, this code is very transparent. It triggered me to create a more compact version:

def imod(a, n):
 c = 1
 while (c % a > 0):
     c += n
 return c // a
Hapten answered 1/4, 2018 at 7:13 Comment(1)
This is OK to explain it to kids, and when n == 7. But otherwise it's about equivalent of this "algorithm": for i in range(2, n): if i * a % n == 1: return iLagena
B
2

from the cpython implementation source code:

def invmod(a, n):
    b, c = 1, 0
    while n:
        q, r = divmod(a, n)
        a, b, c, n = n, c, b - q*c, r
    # at this point a is the gcd of the original inputs
    if a == 1:
        return b
    raise ValueError("Not invertible")

according to the comment above this code, it can return small negative values, so you could potentially check if negative and add n when negative before returning b.

Byronbyrum answered 15/8, 2019 at 1:34 Comment(2)
"so you could potentially check if negative and add n when negative before returning b". Unfortunately n is 0 at that point. (You'd have to save, and use, the original value of n.)Variate
@DonHatch agreedJonquil
N
1

To figure out the modular multiplicative inverse I recommend using the Extended Euclidean Algorithm like this:

def multiplicative_inverse(a, b):
    origA = a
    X = 0
    prevX = 1
    Y = 1
    prevY = 0
    while b != 0:
        temp = b
        quotient = a/b
        b = a%b
        a = temp
        temp = X
        a = prevX - quotient * X
        prevX = temp
        temp = Y
        Y = prevY - quotient * Y
        prevY = temp

    return origA + prevY
Norine answered 21/1, 2012 at 20:33 Comment(1)
There appears to be a bug in this code: a = prevX - quotient * X should be X = prevX - quotient * X, and it should return prevX. FWIW, this implementation is similar to that in Qaz's link in the comment to Märt Bakhoff's answer.Birdhouse
P
0

Well, here's a function in C which you can easily convert to python. In the below c function extended euclidian algorithm is used to calculate inverse mod.

int imod(int a,int n){
int c,i=1;
while(1){
    c = n * i + 1;
    if(c%a==0){
        c = c/a;
        break;
    }
    i++;
}
return c;}

Translates to Python Function

def imod(a,n):
  i=1
  while True:
    c = n * i + 1;
    if(c%a==0):
      c = c/a
      break;
    i = i+1
  
  return c

Reference to the above C function is taken from the following link C program to find Modular Multiplicative Inverse of two Relatively Prime Numbers

Palmieri answered 29/1, 2018 at 14:48 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.