k = 31
F = GF(2)
E = GF(2^k)
PR = PolynomialRing(F, 'x')
p = PR.irreducible_element(k)

while(True):
    gamma = E.random_element()
    if gamma not in F:
        break

# gamma = PR('x^3+x')
# print('0b'+ ''.join(map(str, gamma.coefficients(sparse=False)[::-1])))

eqs = [PR(gamma)^i % p for i in range(k+1)]
for eq in eqs:
    print(eq)
print()

eqs = list(map(lambda eq: eq.coefficients(sparse=False) + [0]*(k-eq.degree()-1), eqs))
for eq in eqs:
    print(eq)
print()
eqs = matrix(eqs[::-1])
print(eqs)
print()
eqs = eqs.transpose()
print(eqs)
print()

rref = eqs.rref()
print(rref)
irrpol = PR(rref[:,-1].list() + [1])
print(irrpol)
print(irrpol.is_irreducible())
print()

# manual rref
def calc_rref_Z2(A, m, n):
    """Calculate the reduced row echelon form of a (mxn)-matrix"""
    A = copy(A)

    # calculate row echelon form
    # https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode
    h = 0 # Initialization of the pivot row
    k = 0 # Initialization of the pivot column

    while h < m and k < n:
        # Find the k-th pivot:
        i_max = -1
        for i in range(h, m):
            if A[i,k] == 1:
                i_max = i
                break

        if i_max == -1:
            # No pivot in this column, pass to next column
            k += 1
        else:
            tmp = A[i_max, :]
            A[i_max, :] = A[h, :]
            A[h, :] = tmp


            # Do for all rows below pivot:
            for i in range(h+1, m):
                f = A[i][k];
                if f == 0:
                    continue;
                A[i, k] = 0;
                for j in range(k+1, n):
                    A[i, j] = A[i, j] - A[h, j]

            # Increase pivot row and column
            h += 1
            k += 1

    # perform back substitution
    for i in range(1, m):
        for j in range(i):
            if A[j, i] == 1:
                A[j, :] -= A[i, :]

    return A

print('Manual RREF:\n')

rref = calc_rref_Z2(eqs, k, k+1) # cols = rows+1, since we transposed the matrix
print(rref)
print()
irrpol = PR(rref[:,-1].list() + [1])
print(irrpol)
print(irrpol.is_irreducible())
print()

print('0b'+ ''.join(map(str, irrpol.coefficients(sparse=False)[::-1])))