La forma más fácil de realizar la inversión de matriz modular con Python?

inverse matrix algorithm c++

Question

Me gustaría tomar la inversa modular de una matriz como [[1,2], [3,4]] mod 7 en Python. He observado numpy (que hace inversión de matriz pero no inversión de matriz modular) y vi algunos paquetes de teoría de números en línea, pero nada que parezca hacer este procedimiento relativamente común (al menos, me parece relativamente común).

Por cierto, el inverso de la matriz anterior es [[5,1], [5,3]] (mod 7). Aunque me gustaría que Python lo haga por mí.

Answer 1

Este pequeño fragmento de código parece hacerlo: enlace

Tenga en cuenta el comentario a continuación para una pequeña mejora. Parece que hace el álgebra lineal correcta hasta donde puedo ver. Nunca he encontrado ninguna opción en los paquetes regulares, por lo que probablemente tomar un fragmento de código de la web (hay muchos más disponibles) es el enfoque más fácil.

Answer 2

De acuerdo ... para los que les importa, resolví mi propio problema. Me tomó un tiempo, pero creo que esto funciona. Probablemente no sea el más elegante, y debería incluir un poco más de manejo de errores, pero funciona:

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 3

Desafortunadamente numpy no tiene implementaciones aritméticas modulares. Siempre puede codificar el algoritmo propuesto usando la reducción de fila o los determinantes como se demuestra aquí . Un inverso modular parece ser bastante útil para la criptografía.

Answer 4

Se puede calcular usando Sage ( www.sagemath.org ) como

Matriz (IntegerModRing (7), [[1, 2], [3,4]]) inversa ()

Aunque Sage es enorme de instalar y tienes que usar la versión de python que viene con él, que es un problema.

Answer 5

Un truco de hack que funciona cuando los errores de redondeo no son un problema:

• encontrar el inverso regular (puede tener entradas no enteras) y el determinante (un entero), ambos implementados en numpy
• multiplica el inverso por el determinante y redondea a enteros (hacky)
• ahora multiplica todo por la inversa multiplicativa del determinante (módulo tu módulo, código debajo)
• hacer mod de entrada por su módulo

Una forma menos hackosa es implementar realmente la eliminación gaussiana. Aquí está mi código usando la eliminación de Gauss, que escribí para mis propios fines (los errores de redondeo fueron un problema para mí). q es el módulo, que no es necesariamente primo.

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

EDITAR: como señalan los comentaristas, hay algunos casos en que este algoritmo falla. No es trivial repararlo, y no tengo tiempo hoy en día. En aquel entonces funcionaba para matrices aleatorias en mi caso (los módulos eran productos de primos grandes). Básicamente, la primera entrada distinta de cero podría no ser relativamente primo para el módulo. El caso principal es fácil ya que puede buscar una fila diferente y cambiar. En el caso no principal, creo que podría ser que todas las entradas principales no son relativamente primarias, por lo que debe combinarlas