Le moyen le plus simple d'effectuer une inversion de matrice modulaire avec Python?

Question

Je voudrais prendre l'inverse modulaire d'une matrice comme [[1,2], [3,4]] mod 7 en Python. J'ai regardé numpy (qui fait l'inversion de matrice mais pas l'inversion de matrice modulaire) et j'ai vu quelques paquets de théorie de nombre en ligne, mais rien qui semble faire cette procédure relativement commune (au moins, il me semble relativement commun). Par ailleurs, l'inverse de la matrice ci-dessus est [[5, 5], [5, 3]] (mod 7). J'aimerais que Python le fasse pour moi.

Answer 1

Ce petit morceau de code semble faire: link

Notez le commentaire ci-dessous pour un peu d'amélioration. Semble faire l'algèbre linéaire correcte autant que je peux voir. Je n'ai jamais trouvé d'option dans les paquets réguliers, donc prendre un extrait de code du web (il y en a beaucoup plus disponible) est l'approche la plus simple.

Answer 2

Malheureusement, numpy ne dispose pas d'implémentations arithmétiques modulaires. Vous pouvez toujours coder l'algorithme proposé en utilisant la réduction de ligne ou les déterminants comme démontré here. Un inverse modulaire semble être très utile pour la cryptographie.

Answer 3

Ok ... pour ceux qui s'en soucient, j'ai résolu mon propre problème. Cela m'a pris du temps, mais je pense que cela fonctionne. Il est sans doute pas le plus élégant, et devrait inclure un peu plus la gestion des erreurs, mais cela fonctionne:

import numpy 
import math 
from numpy import matrix 
from numpy import linalg 

def modMatInv(A,p):  # Finds the inverse of matrix A mod p 
    n=len(A) 
    A=matrix(A) 
    adj=numpy.zeros(shape=(n,n)) 
    for i in range(0,n): 
    for j in range(0,n): 
     adj[i][j]=((-1)**(i+j)*int(round(linalg.det(minor(A,j,i)))))%p 
    return (modInv(int(round(linalg.det(A))),p)*adj)%p 

def modInv(a,p):   # Finds the inverse of a mod p, if it exists 
    for i in range(1,p): 
    if (i*a)%p==1: 
     return i 
    raise ValueError(str(a)+" has no inverse mod "+str(p)) 

def minor(A,i,j): # Return matrix A with the ith row and jth column deleted 
    A=numpy.array(A) 
    minor=numpy.zeros(shape=(len(A)-1,len(A)-1)) 
    p=0 
    for s in range(0,len(minor)): 
    if p==i: 
     p=p+1 
    q=0 
    for t in range(0,len(minor)): 
     if q==j: 
     q=q+1 
     minor[s][t]=A[p][q] 
     q=q+1 
    p=p+1 
    return minor 

Answer 4

Une astuce hackish qui fonctionne lorsque des erreurs d'arrondi ne sont pas un problème:

• trouver l'inverse régulière (peuvent avoir des entrées non-entier), et le déterminant (un entier), tous deux mis en oeuvre dans numpy
• multiplier l'inverse par le déterminant, et autour de nombres entiers (hacky)
• multiplier maintenant tout en inverse multiplicatif du déterminant (modulo votre modulu s, le code ci-dessous)
• n'entrywise mod par votre module

Une manière moins hackish est de mettre effectivement en œuvre l'élimination gaussienne. Voici mon code utilisant l'élimination gaussienne, que j'ai écrit pour mes propres fins (les erreurs d'arrondi étaient un problème pour moi). q est le module, qui n'est pas nécessairement premier.

def generalizedEuclidianAlgorithm(a, b): 
    if b > a: 
     return generalizedEuclidianAlgorithm(b,a); 
    elif b == 0: 
     return (1, 0); 
    else: 
     (x, y) = generalizedEuclidianAlgorithm(b, a % b); 
     return (y, x - (a/b) * y) 

def inversemodp(a, p): 
    a = a % p 
    if (a == 0): 
     print "a is 0 mod p" 
     return None 
    if a > 1 and p % a == 0: 
     return None 
    (x,y) = generalizedEuclidianAlgorithm(p, a % p); 
    inv = y % p 
    assert (inv * a) % p == 1 
    return inv 

def identitymatrix(n): 
    return [[long(x == y) for x in range(0, n)] for y in range(0, n)] 

def inversematrix(matrix, q): 
    n = len(matrix) 
    A = np.matrix([[ matrix[j, i] for i in range(0,n)] for j in range(0, n)], dtype = long) 
    Ainv = np.matrix(identitymatrix(n), dtype = long) 
    for i in range(0, n): 
     factor = inversemodp(A[i,i], q) 
     if factor is None: 
      raise ValueError("TODO: deal with this case") 
     A[i] = A[i] * factor % q 
     Ainv[i] = Ainv[i] * factor % q 
     for j in range(0, n): 
      if (i != j): 
       factor = A[j, i] 
       A[j] = (A[j] - factor * A[i]) % q 
       Ainv[j] = (Ainv[j] - factor * Ainv[i]) % q 
    return Ainv 

EDIT: comme le soulignent les commentateurs, cet algorithme échoue dans certains cas. C'est un peu trivial à corriger, et je n'ai pas le temps de nos jours. À l'époque, cela fonctionnait pour des matrices aléatoires dans mon cas (les modules étaient des produits de grands nombres premiers). Fondamentalement, la première entrée non nulle pourrait ne pas être relativement primaire au module. Le cas principal est facile puisque vous pouvez rechercher une ligne différente et échanger. Dans le cas de non-prime, je pense qu'il pourrait être que toutes entrées principales ne sont pas premiers entre eux, donc vous devez les combiner

Answer 5

Il peut être calculé en utilisant Sage (www.sagemath.org) comme

Matrice (IntegerModRing (7), [[1, 2], [3,4]]).inverse()

Bien que Sage est énorme à installer et vous devez utiliser la version de python qui vient avec ce qui est une douleur.