#Demo modeling example

q = 2
n = 4
K = GF(q)

C1 = random_matrix(K, n, n)
while C1.rank() < n:
    C1 = random_matrix(K, n, n)
C2 = random_matrix(K, n, n)
while C2.rank() < n:
    C2 = random_matrix(K, n, n)

A = GL(n, K).random_element()
B = GL(n, K).random_element()
D1 = A*C1*B
D2 = A*C2*B

R = BooleanPolynomialRing(2*n^2, "x", order='degrevlex')
A_ = matrix(R, n, n, [R.gens()[i] for i in range(n^2)])
B_ = matrix(R, n, n, [R.gens()[i+n^2] for i in range(n^2)])
A_*B_

D1_var = A_*C1*B_
D2_var = A_*C2*B_

system = []
for i in range(n):
    for j in range(n):
        system = system + [D1_var[i,j]] + [D2_var[i,j]]
system


#I = ideal(system)
#B = I.groebner_basis()
#solution = I.variety()

#print(B)
#print(solution)
system