Хакерский трюк, который работает, когда ошибки округления не являются проблемой:
- находит регулярное обратное (может иметь нецелые записи) и определитель (целое число), оба реализованные вnumpy
- умножить обратное на определитель и округлить до целых (хаки)
- теперь умножить все на умножительное обратное определителя (по модулю ваш модуль, код ниже)
- doentrywise mod по вашему модулю
Менее хакерский способ - это реализовать гауссово исключение.Вот мой код с использованием исключения Гаусса, который я написал для своих собственных целей (ошибки округления были проблемой для меня).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
РЕДАКТИРОВАТЬ: как указывают комментаторы, в некоторых случаях этот алгоритм не работает.Это немного нетривиально исправить, и у меня нет времени в настоящее время.Тогда в моем случае это работало для случайных матриц (модули были произведениями больших простых чисел).По сути, первая ненулевая запись не может быть относительно простой по отношению к модулю.Основной случай прост, так как вы можете искать другую строку и менять местами.В не простом случае, я думаю, что все ведущие записи не являются относительно простыми, поэтому вы должны объединить их