Python 中的模乘逆函数-IT科技

摘要：问题描述：某些标准 Python 模块是否包含计算某个数的模乘逆元的函数，例如y = invmod(x, p)满足的数x*y == 1 (mod p)？谷歌似乎没有给出任何好的提示。当然，人们可以想出自制的 10 行扩展欧几里得算法，但为什么要重新发明轮子呢？比如 Java 的BigIntegerhasmo...

问题描述：

某些标准 Python 模块是否包含计算某个数的模乘逆元的函数，例如y = invmod(x, p)满足的数x*y == 1 (mod p)？谷歌似乎没有给出任何好的提示。

当然，人们可以想出自制的 10 行扩展欧几里得算法，但为什么要重新发明轮子呢？

比如 Java 的BigIntegerhasmodInverse方法，Python 不是也有类似的东西吗？

解决方案 1：

Python 3.8+

y = pow(x, -1, p)

Python 3.7 及更早版本

也许有人会发现这很有用（来自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

解决方案 2：

如果模数是素数（你称之为p），那么你可以简单地计算：

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

或者在 Python 中：

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

这里有人用 Python 实现了一些数论功能：http://www.math.umbc.edu/~campbell/Computers/Python/numbthy.html

这是在提示符下完成的一个例子：

m = 1000000007
x = 1234567
y = pow(x,m-2,m)
y
989145189L
x*y
1221166008548163L
x*y % m
1L

解决方案 3：

您可能还想看看gmpy模块。它是 Python 和 GMP 多精度库之间的接口。gmpy 提供了一个 invert 函数，可以完全满足您的需求：

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

更新答案

正如@hyh所指出的，gmpy.invert()如果逆不存在，则返回0。这与GMPmpz_invert()函数的行为一致。gmpy.divm(a, b, m)提供了的通用解决方案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()当时将返回一个解决方案gcd(b,m) == 1，当乘法逆元不存在时将引发异常。

免责声明：我是 gmpy 库的当前维护者。

更新答案 2

当逆不存在时，gmpy2 现在会正确地引发异常：

>>> import gmpy2

>>> gmpy2.invert(0,5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: invert() no inverse exists

解决方案 4：

从 3.8 开始，python 的 pow() 函数可以接受模数和负整数。参见此处。其使用方法如下：

>>> pow(38, -1, 97)
23
>>> 23 * 38 % 97 == 1
True

解决方案 5：

这是CodeFights的一行代码；它是最短的解决方案之一：

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]

-1如果A没有乘法逆元，它将返回n。

用法：

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

该解决方案采用扩展欧几里得算法。

解决方案 6：

Sympy是一个用于符号数学的 Python 模块，如果您不想实现自己的模块逆函数（或者您已经在使用 Sympy），它有一个内置的模块逆函数：

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'

Sympy 网站上似乎没有记录这一点，但这里有文档字符串：Github 上的 Sympy mod_inverse 文档字符串

解决方案 7：

这是一个简洁的单行代码，无需使用任何外部库即可完成此操作。

# 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

请注意，这实际上只是 egcd，简化为仅返回感兴趣的单个系数。

解决方案 8：

我尝试了该线程中的不同解决方案，最后我使用了这个：

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

Python中的Modular_inverse

解决方案 9：

这是我的代码，它可能有点草率，但无论如何它似乎对我有用。

# 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

解决方案 10：

来自 cpython 实现源代码：

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")

根据此代码上方的注释，它可以返回较小的负值，因此您可以在返回 b 之前检查是否为负数，并在为负数时添加 n。

解决方案 11：

为了计算模乘逆元，我建议使用扩展欧几里得算法，如下所示：

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

解决方案 12：

上面的代码无法在 Python 3 中运行，而且与 GCD 版本相比效率较低。不过，这段代码非常清晰易懂。这促使我创建了一个更紧凑的版本：

def imod(a, n):
 c = 1
 while (c % a > 0):
     c += n
 return c // a

解决方案 13：

好吧，这里有一个 C 函数，你可以轻松地将其转换为 Python 代码。下面的 C 函数使用了扩展欧几里得算法来计算逆模。

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;}

转换为 Python 函数

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

上述 C 函数的引用取自以下链接：C 程序用于查找两个相对素数的模乘逆元