#!/usr/bin/env python3 -I
"""
Search (char 2) for f = x^2 + eps(y,z,w), eps a sum of monomials, and a monomial
curve C : t -> (t^A, t^a, t^b, t^c) in A^4 such that

  (1) f vanishes on C,
  (2) d_y eps, d_z eps, d_w eps vanish on C          (so ord_eta f >= 2 on C),
  (3) the mixed Hasse derivatives d_y d_z, d_y d_w, d_z d_w of eps vanish on C
      (so in_eta(f) = (x - x_eta)^2 + sum_j c_j (y_j - eta_j)^2 is the square of a
       linear form at every eta in C: the "edge data" of the manuscript is
       (n, n-r, q) = (4, 3, 2) at every point of C),
  (4) the semigroup <A, a, b, c> has all four of A, a, b, c as minimal generators
      (so C has embedding dimension 4 at the origin and is not contained in
       any smooth hypersurface through 0),
  (5) every monomial of eps has total degree >= 3 (so in_0 f = x^2).

Parity bookkeeping (char 2): for monomial m = y^i z^j w^k of t-weight
D(m) = i a + j b + k c and parity class (i,j,k) mod 2:
  f|_C = 0            <=>  for each D, #{m : D(m)=D} = [D == 2A]   (mod 2)
  d_y eps|_C = 0      <=>  for each D, #{m : D(m)=D, i odd} = 0     (mod 2)
  d_y d_z eps|_C = 0  <=>  for each D, #{m : D(m)=D, i,j odd} = 0   (mod 2)
  etc.  Solving the 7 linear conditions on the 8 class counts (mod 2) gives:
  at weight D != 2A all eight class counts have the same parity; at D = 2A the
  seven nonzero classes are odd and class 000 is even (the alternative,
  class 000 odd and the rest even, forces C into the smooth hypersurface
  x = monomial).  Minimal solution: one monomial in each of the 7 nonzero
  parity classes, all of weight 2A.
"""
import itertools, sys

def minimal_generators(gens):
    gens = sorted(set(gens))
    # semigroup membership up to max
    M = 2 * max(gens) + 1
    reach = [False] * (M + 1)
    reach[0] = True
    for n in range(1, M + 1):
        for g in gens:
            if n - g >= 0 and reach[n - g]:
                reach[n] = True
                break
    mins = []
    for g in gens:
        others = [h for h in gens if h != g]
        # is g in semigroup generated by others?
        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 search(maxw=40, maxdeg=12):
    found = []
    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:
                        continue
                    import math
                    if math.gcd(math.gcd(A, a), math.gcd(b, c)) != 1:
                        continue
                    if len(minimal_generators(gens)) < 4:
                        continue
                    target = 2 * A
                    # monomials of weight 2A in each nonzero parity class, degree >= 3
                    classes = {cl: [] for cl in itertools.product((0, 1), repeat=3) if cl != (0, 0, 0)}
                    for i in range(0, maxdeg + 1):
                        for j in range(0, maxdeg + 1 - i):
                            for k in range(0, maxdeg + 1 - i - j):
                                if i + j + k < 3:
                                    continue
                                if i * a + j * b + k * c == target:
                                    cl = (i % 2, j % 2, k % 2)
                                    if cl != (0, 0, 0):
                                        classes[cl].append((i, j, k))
                    if all(classes[cl] for cl in classes):
                        found.append((A, a, b, c, {cl: classes[cl] for cl in classes}))
                        if len(found) >= 40:
                            return found
    return found

if __name__ == "__main__":
    res = search()
    print("number of (A,a,b,c) with a 7-monomial solution:", len(res))
    for A, a, b, c, classes in res[:40]:
        print(f"\nC = (t^{A}, t^{a}, t^{b}, t^{c}),  2A = {2*A}")
        for cl, ms in sorted(classes.items()):
            print("   class", cl, "->", ms[:4])