Fungsi invers perkalian modular dengan Python

110

Apakah beberapa modul Python standar berisi fungsi untuk menghitung invers perkalian modular dari sebuah bilangan, yaitu bilangan y = invmod(x, p)seperti itu x*y == 1 (mod p)? Google sepertinya tidak memberikan petunjuk bagus tentang ini.

Tentu saja, seseorang dapat membuat 10-liner buatan sendiri dari algoritma Euclidean yang diperluas , tetapi mengapa menemukan kembali roda.

Misalnya, Java BigIntegermemiliki modInversemetode. Bukankah Python memiliki sesuatu yang serupa?

dorserg
sumber
18
Python 3.8 (yang akan dirilis akhir tahun ini), Anda akan dapat menggunakan built-in powfungsi untuk ini: y = pow(x, -1, p). Lihat bugs.python.org/issue36027 . Hanya butuh 8,5 tahun dari pertanyaan yang diajukan hingga solusi muncul di perpustakaan standar!
Mark Dickinson
4
Saya melihat @MarkDickinson dengan rendah hati tidak menyebutkan bahwa ey adalah pembuat peningkatan yang sangat berguna ini, jadi saya akan melakukannya. Terima kasih untuk pekerjaan ini, Mark, ini tampak hebat!
Don Hatch

Jawaban:

128

Mungkin seseorang akan menganggap ini berguna (dari wikibook ):

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
Märt Bakhoff
sumber
1
Saya mengalami masalah dengan bilangan negatif menggunakan algoritma ini. modinv (-3, 11) tidak berfungsi. Saya memperbaikinya dengan mengganti egcd dengan penerapan di halaman dua pdf ini: anh.cs.luc.edu/331/notes/xgcd.pdf Semoga membantu!
Qaz
@Qaz Anda juga dapat mengurangi -3 modulo 11 untuk membuatnya positif, dalam hal ini modinv (-3, 11) == modinv (-3 + 11, 11) == modinv (8, 11). Mungkin itulah yang kebetulan dilakukan algoritme dalam PDF Anda di beberapa titik.
Thomas
1
Jika Anda kebetulan menggunakan sympy, x, _, g = sympy.numbers.igcdex(a, m)lakukan triknya.
Lynn
59

Jika modulus Anda adalah bilangan prima (Anda menyebutnya p) maka Anda dapat menghitung:

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

Atau dengan Python:

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

Berikut adalah seseorang yang telah menerapkan beberapa kemampuan teori bilangan dengan Python: http://www.math.umbc.edu/~campbell/Computers/Python/numbthy.html

Berikut adalah contoh yang dilakukan saat diminta:

m = 1000000007
x = 1234567
y = pow(x,m-2,m)
y
989145189L
x*y
1221166008548163L
x*y % m
1L
phkahler.dll
sumber
1
Eksponen yang naif bukanlah pilihan karena batas waktu (dan memori) untuk nilai p yang cukup besar seperti misalnya 1000000007.
dorserg
16
eksponensial modular dilakukan paling banyak dengan perkalian N * 2 di mana N adalah jumlah bit dalam eksponen. menggunakan modulus 2 ** 63-1 invers dapat dihitung pada prompt dan segera mengembalikan hasilnya.
phkahler
3
Wow Keren. Saya menyadari eksponensial cepat, saya hanya tidak menyadari bahwa fungsi pow () dapat mengambil argumen ketiga yang mengubahnya menjadi eksponensial modular.
dorserg
5
Itu sebabnya Anda menggunakan Python, kan? Karena itu luar biasa :-)
phkahler
2
Ngomong-ngomong ini bekerja karena dari teorema Fermat pow kecil (x, m-1, m) harus 1. Jadi (pow (x, m-2, m) * x)% m == 1. Jadi pow (x, m-2, m) adalah kebalikan dari x (mod m).
Piotr Dabkowski
21

Anda mungkin juga ingin melihat modul gmpy . Ini adalah antarmuka antara Python dan pustaka multi-presisi GMP. gmpy menyediakan fungsi invert yang melakukan apa yang Anda butuhkan:

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

Jawaban yang diperbarui

Seperti dicatat oleh @hyh, gmpy.invert()mengembalikan 0 jika invers tidak ada. Itu cocok dengan perilaku mpz_invert()fungsi GMP . gmpy.divm(a, b, m)memberikan solusi umum untuk 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()akan mengembalikan solusi saat gcd(b,m) == 1dan memunculkan pengecualian saat pembalikan perkalian tidak ada.

Penafian: Saya adalah pengelola pustaka gmpy saat ini.

Jawaban yang diperbarui 2

gmpy2 sekarang memunculkan pengecualian dengan benar saat inversi tidak ada:

>>> import gmpy2

>>> gmpy2.invert(0,5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ZeroDivisionError: invert() no inverse exists
casevh
sumber
Ini keren sampai saya menemukan gmpy.invert(0,5) = mpz(0)alih-alih memunculkan kesalahan ...
h__
@hyh Bisakah Anda melaporkan ini sebagai masalah di beranda gmpy? Itu selalu dihargai jika masalah dilaporkan.
casevh
BTW, apakah ada perkalian modular dalam gmpypaket ini ? (yaitu beberapa fungsi yang memiliki nilai yang sama tetapi lebih cepat dari (a * b)% p?)
h__
Ini telah diusulkan sebelumnya dan saya sedang bereksperimen dengan metode yang berbeda. Pendekatan paling sederhana untuk hanya menghitung (a * b) % pdalam suatu fungsi tidak lebih cepat dari sekedar mengevaluasi (a * b) % pdengan Python. Overhead untuk panggilan fungsi lebih besar daripada biaya evaluasi ekspresi. Lihat code.google.com/p/gmpy/issues/detail?id=61 untuk detail selengkapnya.
casevh
2
Hal yang hebat adalah bahwa ini juga berfungsi untuk modulus non-prima.
synecdoche
13

Pada 3.8 pythons pow () fungsi dapat mengambil modulus dan integer negatif. Lihat disini . Kasus mereka tentang cara menggunakannya adalah

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

Berikut adalah satu baris untuk CodeFights ; ini adalah salah satu solusi terpendek:

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]

Ini akan kembali -1jika Atidak memiliki pembalikan perkaliann .

Pemakaian:

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

Solusinya menggunakan Algoritma Euclidean yang Diperluas .

HKTonyLee
sumber
6

Sympy , modul python untuk matematika simbolik, memiliki fungsi invers modular bawaan jika Anda tidak ingin mengimplementasikannya sendiri (atau jika Anda sudah menggunakan 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'

Ini sepertinya tidak didokumentasikan di situs web Sympy, tetapi inilah docstringnya: Sympy mod_inverse docstring di Github

Chris Chudzicki
sumber
2

Ini kode saya, mungkin ceroboh tetapi tampaknya tetap berfungsi untuk saya.

# 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
Eric
sumber
2

Kode di atas tidak akan berjalan di python3 dan kurang efisien dibandingkan dengan varian GCD. Namun, kode ini sangat transparan. Ini memicu saya untuk membuat versi yang lebih ringkas:

def imod(a, n):
 c = 1
 while (c % a > 0):
     c += n
 return c // a
BvdM
sumber
1
Tidak apa-apa menjelaskannya kepada anak-anak, dan kapan n == 7. Tetapi sebaliknya, ini hampir setara dengan "algoritme" ini:for i in range(2, n): if i * a % n == 1: return i
Tomasz Gandor
2

Berikut adalah 1-liner ringkas yang melakukannya, tanpa menggunakan library eksternal apa pun.

# 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

Perhatikan bahwa ini benar-benar hanya egcd, disederhanakan untuk mengembalikan hanya koefisien tunggal yang menarik.

Don Hatch
sumber
1

Untuk mengetahui pembalikan perkalian modular saya sarankan menggunakan Algoritma Euclidean yang Diperluas seperti ini:

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
David Sulpy
sumber
Tampaknya ada bug dalam kode ini: a = prevX - quotient * Xseharusnya X = prevX - quotient * X, dan seharusnya kembali prevX. FWIW, implementasi ini mirip dengan yang ada di link Qaz di komentar jawaban Märt Bakhoff.
PM 2Ring
1

Saya mencoba solusi yang berbeda dari utas ini dan pada akhirnya saya menggunakan yang ini:

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 dengan Python

Al Po
sumber
1
kode ini tidak valid. returndi egcd diindendasikan dengan cara yang salah
ph4r05
0

Yah, saya tidak memiliki fungsi di python tetapi saya memiliki fungsi di C yang dapat Anda ubah dengan mudah ke python, di bawah fungsi c, algoritma euclidian diperpanjang digunakan untuk menghitung mod terbalik.

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

Fungsi 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

Referensi fungsi C di atas diambil dari link program C berikut untuk mencari Pembalikan Perkalian Modular dari dua Bilangan Prima Relatif

Mohd Shibli
sumber
0

dari kode sumber implementasi 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")

menurut komentar di atas kode ini, ia dapat mengembalikan nilai negatif kecil, jadi Anda berpotensi memeriksa apakah negatif dan menambahkan n ketika negatif sebelum mengembalikan b.

micsthepick
sumber
"sehingga Anda berpotensi memeriksa apakah negatif dan menambahkan n ketika negatif sebelum mengembalikan b". Sayangnya n adalah 0 pada saat itu. (Anda harus menyimpan, dan menggunakan, nilai asli n.)
Don Hatch