python求梯形面积,Python：非常大的矩阵的简化行梯形形式（mod p）

I want to find find a reduced a row echelon form (in field F_q) of a big matrix.

I tried the following code.

Although I used gmpy2 library to speed up, the program was still out of memory. because my input matrix is very large (100 x 2^15) and p is also very large (|p|=256 bits). Can someone suggest how to reduce the complexity of this alg.

Thank you

def invmodp(a, p):

return gmpy2.invert(a,p)

def division_mod(a, b, p): #a/b mod p

invert = invmodp(b, p)

return (a * invert) %p

def row_echelon_form(M, p):

lead = 0

rowCount = len(M)

columnCount = len(M[0])

for r in range(rowCount):

if lead >= columnCount:

return

i = r

while M[i][lead] == 0:

i += 1

if i == rowCount:

i = r

lead += 1

if columnCount == lead:

return

M[i],M[r] = M[r],M[i]

lv = M[r][lead]

M[r] = [ division_mod(mrx, lv, p) for mrx in M[r]]

for i in range(rowCount):

if i != r:

lv = M[i][lead]

M[i] = [ (iv - lv*rv)%p for rv,iv in zip(M[r],M[i])]

lead += 1

return M

解决方案

I was able to save a few seconds of running time by using gmpy2.divm to replace your division_mod. I wasn't able to make any other significant improvements. The following program creates a random 100 x 2^15 matrix and calculates the row echelon form in approximately 3 minutes and consumes 425MB of memory.

import gmpy2

bits = 256

r = 100

c = 2**15

p = gmpy2.next_prime(2**bits - 1234)

seed = gmpy2.random_state(42)

M = []

for i in range(r):

M.append([gmpy2.mpz_urandomb(seed, bits) for j in range(c)])

def row_echelon_form(M, p):

lead = 0

rowCount = len(M)

columnCount = len(M[0])

for r in range(rowCount):

if lead >= columnCount:

return

i = r

while M[i][lead] == 0:

i += 1

if i == rowCount:

i = r

lead += 1

if columnCount == lead:

return

M[i],M[r] = M[r],M[i]

lv = M[r][lead]

M[r] = [ gmpy2.divm(mrx, lv, p) for mrx in M[r]]

for i in range(rowCount):

if i != r:

lv = M[i][lead]

M[i] = [ (iv - lv*rv) % p for rv,iv in zip(M[r],M[i])]

lead += 1

return M

N = row_echelon_form(M, p)

If your memory usage grows beyond about 500MB, there may be a memory leak in your version of gmpy2. Or I've interpreted your requirements incorrectly and the matrix is significantly larger.

Disclaimer: I maintain gmpy2.