from sympy import symbols, Eq, GroebnerBasis

Define variables

x = symbols('x_AA x_AG x_AC x_AT x_GA x_GG x_GC x_GT x_CA x_CG x_CC x_CT x_TA x_TG x_TC x_TT')

Define polynomial equations

equations = [ Eq(x[1] - x[4], 0), # x_AG - x_GA = 0 Eq(x[2] - x[8], 0), # x_AC - x_CA = 0 Eq(x[3] - x[12], 0), # x_AT - x_TA = 0 Eq(x[6] - x[9], 0), # x_GC - x_CG = 0 Eq(x[7] - x[13], 0), # x_GT - x_TG = 0 Eq(x[10] - x[14], 0), # x_CT - x_TC = 0 Eq(x[0] + x[1] + x[2] + x[3] - 1, 0), # x_AA + x_AG + x_AC + x_AT = 1 Eq(x[4] + x[5] + x[6] + x[7] - 1, 0), # x_GA + x_GG + x_GC + x_GT = 1 Eq(x[8] + x[9] + x[10] + x[11] - 1, 0), # x_CA + x_CG + x_CC + x_CT = 1 Eq(x[12] + x[13] + x[14] + x[15] - 1, 0), # x_TA + x_TG + x_TC + x_TT = 1 Eq(x[0] - x[5], 0), # x_AA - x_GG = 0 Eq(x[2] - x[7], 0), # x_AC - x_GT = 0 Eq(x[8] - x[13], 0), # x_CA - x_TG = 0 ]

Generate the ideal using Groebner basis

basis = GroebnerBasis([eq.lhs for eq in equations], x)

Display the result

print(basis) give me the code for generators