#!/usr/bin/env python3 -I
"""
Full F_2-linear search for f = x^2 + eps(y,z,w) in char 2 with a monomial curve
C = (t^A, t^a, t^b, t^c) of embedding dimension 4 such that
  f|_C = 0, (d_y eps, d_z eps, d_w eps)|_C = 0, (d_y d_z, d_y d_w, d_z d_w eps)|_C = 0,
with every monomial of eps of degree >= 3.

Unknowns: coefficients (in F_2) of monomials y^i z^j w^k, 3 <= i+j+k <= maxdeg.
For each t-weight D occurring, 7 equations:
  sum_{m: D(m)=D} c_m                      = [D == 2A]
  sum_{m: D(m)=D, i odd} c_m               = 0
  sum_{m: D(m)=D, j odd} c_m               = 0
  sum_{m: D(m)=D, k odd} c_m               = 0
  sum_{m: D(m)=D, i,j odd} c_m             = 0
  sum_{m: D(m)=D, i,k odd} c_m             = 0
  sum_{m: D(m)=D, j,k odd} c_m             = 0
We solve by Gaussian elimination over F_2 and print a sparse solution.
"""
import itertools, math, sys

def minimal_generators(gens):
    mins = []
    for g in gens:
        others = [h for h in gens if h != g]
        r = [False] * (g + 1); r[0] = True
        for n in range(1, g + 1):
            for h in others:
                if n - h >= 0 and r[n - h]:
                    r[n] = True; break
        if not r[g]:
            mins.append(g)
    return mins

def solve_f2(rows, ncols):
    """rows: list of (bitmask over columns, rhs bit). Returns a solution bitmask or None."""
    pivots = []  # list of (pivot_col, mask, rhs)
    for mask, rhs in rows:
        for pc, pm, pr in pivots:
            if mask >> pc & 1:
                mask ^= pm; rhs ^= pr
        if mask == 0:
            if rhs:
                return None
            continue
        pc = mask.bit_length() - 1
        # reduce existing pivots
        newp = []
        for qc, qm, qr in pivots:
            if qm >> pc & 1:
                qm ^= mask; qr ^= rhs
            newp.append((qc, qm, qr))
        pivots = newp + [(pc, mask, rhs)]
    sol = 0
    for pc, pm, pr in pivots:
        # free variables set to 0: value of pivot var = rhs
        if pr:
            sol |= 1 << pc
    return sol

def try_curve(A, a, b, c, maxdeg):
    monos = []
    for i in range(maxdeg + 1):
        for j in range(maxdeg + 1 - i):
            for k in range(maxdeg + 1 - i - j):
                if i + j + k >= 3:
                    monos.append((i, j, k))
    wt = {m: m[0] * a + m[1] * b + m[2] * c for m in monos}
    idx = {m: n for n, m in enumerate(monos)}
    byw = {}
    for m in monos:
        byw.setdefault(wt[m], []).append(m)
    rows = []
    for D, ms in byw.items():
        pats = [lambda m: True,
                lambda m: m[0] % 2 == 1, lambda m: m[1] % 2 == 1, lambda m: m[2] % 2 == 1,
                lambda m: m[0] % 2 == 1 and m[1] % 2 == 1,
                lambda m: m[0] % 2 == 1 and m[2] % 2 == 1,
                lambda m: m[1] % 2 == 1 and m[2] % 2 == 1]
        for n, pat in enumerate(pats):
            mask = 0
            for m in ms:
                if pat(m):
                    mask |= 1 << idx[m]
            rhs = 1 if (n == 0 and D == 2 * A) else 0
            if mask or rhs:
                rows.append((mask, rhs))
    if 2 * A not in byw:
        return None
    sol = solve_f2(rows, len(monos))
    if sol is None:
        return None
    chosen = [m for m in monos if sol >> idx[m] & 1]
    return chosen

def main(maxw=24, maxdeg=10, nmax=12):
    found = 0
    for A in range(2, maxw):
        for a in range(2, maxw):
            for b in range(a + 1, maxw):
                for c in range(b + 1, maxw):
                    gens = [A, a, b, c]
                    if len(set(gens)) < 4 or math.gcd(math.gcd(A, a), math.gcd(b, c)) != 1:
                        continue
                    if len(minimal_generators(gens)) < 4:
                        continue
                    ch = try_curve(A, a, b, c, maxdeg)
                    if ch:
                        found += 1
                        print(f"C=(t^{A},t^{a},t^{b},t^{c})  eps = " + " + ".join(
                            f"y^{i}z^{j}w^{k}" for i, j, k in ch))
                        if found >= nmax:
                            return
    print("found", found)

if __name__ == "__main__":
    main()