Trying to get the right inverse, not possible

edit

Hey guys,

I am working on a project and I am trying to find a matrix that is invertible using the following code:

e = 5;
n = 2;
m = 3;
s = 2;
t = 2;
b = 2;


mon = (b * n^s)^t
mon2 = b * n^s
#F_q is F_2/{irreducible element in F_2}
F.<r> = GF(2)[];

for p in F.polynomials(e):
    if p.is_irreducible():
        break;

K.<q> = GF(2^e, name='q', modulus=p);

Zn = Integers(2^(e*n)-1);
Zm = Integers(2^(e*m)-1);

R = PolynomialRing(K,'X');
R.inject_variables();

M1 = matrix(K,mon,mon);
M2 = matrix(K,mon,mon);
M1inv = matrix(K,mon,mon);
M2inv = matrix(K,mon,mon);
pt_sec2pub = matrix(K,mon,m*n);
#Constructing matrix A using the variables mentioned in the paper
A = matrix(Zn,m,m);

while True:
    for i in range(0,m):
        for j in range(0,m):
            if(s<m):
                if j==s-i:
                    A[i,j] = 0;
                else:
                    A[i,j] = 2^(ZZ.random_element(0,n*e));
            else:
                A[i,j] = 2^(ZZ.random_element(0,n*e));
    det = A.determinant();
    if gcd(det,2^(e*n)-1) == 1:
        break;
Ainv = matrix(Zn,m,m);
Ainv = A.inverse();

print A * Ainv
print Ainv * A
#Constructing matrix B using the variables mentioned in the paper
B = matrix(Zm,n,n);

while True:
    for i in range(0,n):
        for j in range(0,n):
            if(t<n):
                if j==t-i:
                    B[i,j] = 0;
                else:
                    B[i,j] = 2^(ZZ.random_element(0,m*e));
            else:
                B[i,j] = 2^(ZZ.random_element(0,m*e));
    det = B.determinant();
    if gcd(det,2^(e*m)-1) == 1:
        break;
Binv = matrix(Zm,m,m);
Binv = B.inverse();
print Binv * B
print B * Binv

# This computes the monomials generated when you apply G1 and G2 to the input. Instead of raising polynomials, you raise the elements of the polynomials to their respective element in A and B.
def compute_monomials_for_public_key(x):
    row = matrix(K,t,mon2)
    vector_result = matrix(K,mon*2,1)

    # This is the vec1^(A[0][0]).lift() * vec2^(A[0][1]).lift() from G1
    row[0,0] = x[0][0]^A[0][0].lift() * x[2][0] ^ A[0][1].lift()
    row[0,1] = x[1][0]^A[0][0].lift() * x[2][0] ^ A[0][1].lift()
    row[0,2] = x[0][0]^A[0][0].lift() * x[3][0] ^ A[0][1].lift()
    row[0,3] = x[1][0]^A[0][0].lift() * x[3][0] ^ A[0][1].lift()
    # This is the vec1^(A[1][0]).lift() * vec3^(A[1][2]).lift(); from G1
    row[0,4] = x[0][0]^A[1][0].lift() * x[4][0] ^ A[1][2].lift()
    row[0,5] = x[1][0]^A[1][0].lift() * x[4][0] ^ A[1][2].lift()
    row[0,6] = x[0][0]^A[1][0].lift() * x[5][0] ^ A[1][2].lift()
    row[0,7] = x[1][0]^A[1][0].lift() * x[5][0] ^ A[1][2].lift()

    # This is the vec1^(A[1][0]).lift() * vec3^(A[1][2]).lift(); from G1
    row[1,0] = x[0][0]^A[1][0].lift() * x[4][0] ^ A[1][2].lift()
    row[1,1] = x[1][0]^A[1][0].lift() * x[4][0] ^ A[1][2].lift()
    row[1,2] = x[0][0]^A[1][0].lift() * x[5][0] ^ A[1][2].lift()
    row[1,3] = x[1][0]^A[1][0].lift() * x[5][0] ^ A[1][2].lift()
    # This is the vec2^(A[2][1]).lift() * vec3^(A[2][2]).lift(); from G1
    row[1,4] = x[2][0]^A[2][1].lift() * x[4][0] ^ A[2][2].lift()
    row[1,5] = x[3][0]^A[2][1].lift() * x[4][0] ^ A[2][2].lift()
    row[1,6] = x[2][0]^A[2][1].lift() * x[5][0] ^ A[2][2].lift()
    row[1,7] = x[3][0]^A[2][1].lift() * x[5][0] ^ A[2][2].lift()

    #print row
    for i in range(0,mon2):
        for j in range(0,mon2):
            vector_result[mon2*i + j] = row[0][i] ^ B[0][0].lift() * row[1][j] ^ B[0][1].lift()
            vector_result[mon2*i + j + mon] = row[0][i] ^ B[1][0].lift() * row[1][j] ^ B[1][1].lift()

    return vector_result




def generateVectors():
    while True:
        for i in range(0,mon):
            for j in range(0,m*n):
                pt_sec2pub[i,j] = K.random_element();
            transp = matrix(pt_sec2pub[i]).transpose();
            #print "HA"
            vec_transp = compute_monomials_for_public_key(transp);
            for j in range(0,mon):
                M1[j,i] = vec_transp[j][0];
                M2[j,i] = vec_transp[j+mon][0];
            #if i == 0:
                #print vec_transp
        print M1.is_invertible()
        M1inv = M1.inverse()
        print M1inv * M1
        break;

generateVectors()

However, when e < 8, it always gives me that M1 is not invertible and I do not understand whether there is an issue in my logic or not. I am however able to compute the inverse of M1 but when I multiply it with M1, I expect the identity matrix, however this is not the case when e < 8. Please try with different values of e to see this happening. Please let me know if you manage to find anything.

Thanks