我想在Python中使用[[1,2],[3,4]] mod 7等矩陣的模逆。我已經看過numpy(它進行矩陣求逆但不是模塊矩陣求逆),我在線看到了一些數字理論包,但似乎沒有做這個相對常見的過程(至少對我來說似乎比較常見)。用Python執行模塊化矩陣反演的最簡單方法是什麼?
順便說一下,上述矩陣的逆是[[5,1],[5,3]](mod 7)。我希望Python爲我做到這一點。
我想在Python中使用[[1,2],[3,4]] mod 7等矩陣的模逆。我已經看過numpy(它進行矩陣求逆但不是模塊矩陣求逆),我在線看到了一些數字理論包,但似乎沒有做這個相對常見的過程(至少對我來說似乎比較常見)。用Python執行模塊化矩陣反演的最簡單方法是什麼?
順便說一下,上述矩陣的逆是[[5,1],[5,3]](mod 7)。我希望Python爲我做到這一點。
不幸的是numpy沒有模塊化的算術實現。您可以總是使用行減少或行列式對提議的算法進行編碼,如演示here。模塊化反演似乎對密碼學非常有用。
好的......對於那些關心的人,我解決了我自己的問題。我花了一段時間,但我認爲這是有效的。它可能不是最優雅,而應該包括一些錯誤處理,但它的工作原理:
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
一個hackish的伎倆舍入誤差時的作品是不是一個問題:
一個不太hackish的方法是實際執行高斯消元法國防部。這裏是我使用高斯消元的代碼,這是我爲自己的目的而寫的(舍入誤差對我來說是一個問題)。 q是模數,不一定是素數。
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
編輯:作爲評論者指出,有一些情況下,這種算法失敗。這個問題稍微不重要,現在我沒時間了。那時它在我的情況下適用於隨機矩陣(模數是大質數的產物)。基本上,第一個非零入口可能不是模數的相對質數。主要情況很簡單,因爲您可以搜索不同的行並進行交換。在非黃金的情況下,我認爲這可能是所有領先的條目不互質,所以你必須把它們混合起來
它可以通過賢者(www.sagemath.org)作爲
計算矩陣(IntegerModRing(7),[[1,2],[3,4]])。逆()
雖然聖人是巨大的安裝,你必須使用其附帶這是一種痛苦的Python版本。
你看看賢者嗎? – Alejandro 2010-11-27 03:22:09
如果你最終編寫自己的一小段代碼。請考慮在這裏分享它,因爲我認爲我們很多人可能會感興趣:)。 – bastijn 2010-11-27 12:42:08
模塊化矩陣求逆被嵌入到'sympy`中(可能是新的,因爲這個問題被問到),模塊化的行減少也相當容易。請參閱http://stackoverflow.com/a/37015283/2747370 – 2016-05-05 16:05:10