# Coordinator follow-up for the harvest of
# ai/2026-10-23-hironaka-char-p-resolution-review (SageMath 10.8).
# Family E = (x^p + y^n, p) over GF(p), p prime, p does not divide n, n > 2p.
# Checks, for p in {2,3,5,7} and several n:
#  (a) Sing(E) = {0}: radical of the ideal of Hasse derivatives of order < p is (x,y);
#  (b) the Hasse derivative d_x^(p) f = 1 (so wp^nega(E,-a) = O, Thm A);
#  (c) blowup of 0: chart y gives x^p + y^(n-p), chart x has no point of order >= p
#      over GF(p^2) near the exceptional divisor (brute force on the exceptional line);
#  (d) monomial valuation v(x)=n, v(y)=p: min over generators (Diff^(j) f in degree p-j)
#      of v/degree equals n;
#  (e) arc x = -s^n + s^N, y = s^p (N large): f(arc) = s^(Np); min over generators
#      of w/degree equals min(N, p(n-1)/(p-1)); and n*d < d*p(n-1)/(p-1) for d>=1;
#  (f) AR reduction to {x=0} (Thm 3.10): J#|_{x=0} = (y^(n (p-1)!)), with exponent b# = p!;
#  (g) the same exponent for n-p equals n(p-1)! - p!, and p does not divide either.
import itertools
from math import factorial
def hasse(f, beta):
    R = f.parent(); K = R.base_ring(); res = R(0)
    for mon, c in f.dict().items():
        e = list(mon); coef = K(c); ok = True
        for i in range(R.ngens()):
            if e[i] < beta[i]: ok = False; break
            coef *= K(binomial(e[i], beta[i])); e[i] -= beta[i]
        if ok and coef != 0: res += coef*R.monomial(*e)
    return res
def ordat(f, pt):
    R = f.parent(); g = f(*[R.gen(i)+pt[i] for i in range(R.ngens())])
    return oo if g == 0 else min(sum(m) for m in g.dict())
allok = True
for p in [2,3,5,7]:
    ns = [n for n in range(2*p+1, 2*p+12) if n % p != 0][:4]
    for n in ns:
        R.<x,y> = PolynomialRing(GF(p)); f = x^p + y^n
        betas = [b for b in itertools.product(range(p), repeat=2) if sum(b) < p]
        sing = R.ideal([hasse(f,b) for b in betas]).radical()
        a_ok = (sing == R.ideal([x,y]))
        b_ok = (hasse(f,(p,0)) == 1)
        f1 = R(f(x*y, y) // y^p); c_ok = (f1 == x^p + y^(n-p))
        F2 = GF(p^2,'a'); R2 = R.change_ring(F2)
        g2 = R2(f(x, y*x) // x^p)   # chart x: exceptional divisor x = 0
        c2 = all(ordat(g2, (F2(0), t)) < p for t in F2)
        c1 = all(ordat(R2(f1), (t, F2(0))) < p for t in F2 if t != 0) and ordat(f1,(0,0)) == p
        # generators of A: Diff^(j) f in degree p - j, j = 0..p-1 (Hasse derivatives of order <= j)
        gens = []
        for j in range(p):
            for b in itertools.product(range(j+1), repeat=2):
                if sum(b) <= j:
                    h = hasse(f, b)
                    if h != 0: gens.append((h, p-j))
        def vmono(h):  # monomial valuation v(x)=n, v(y)=p
            return min(n*m[0] + p*m[1] for m in h.dict())
        mu_mono = min(QQ(vmono(h))/d for h,d in gens)
        d_ok = (mu_mono == n)
        N = 2*n + 3   # any N > p(n-1)/(p-1) works
        S.<s> = PowerSeriesRing(GF(p), default_prec=N*p+10)
        arc = (-s^n + s^N, s^p)
        def w(h): return h(*arc).valuation()
        e1 = (w(f) == N*p)
        mu_arc = min(QQ(w(h))/d for h,d in gens)
        e2 = (mu_arc == min(N, QQ(p*(n-1))/(p-1)))
        e3 = all(n*d < d*mu_arc for d in range(1, p))
        # AR reduction
        bs = factorial(p); ej = []
        for j in range(p):
            # restrict Diff^(j) f to x = 0; in K[y] it is principal, generated by y^(min order)
            Dj = [R(hasse(f,b).subs(x=0)) for b in itertools.product(range(j+1), repeat=2) if sum(b) <= j]
            Dj = [h for h in Dj if h != 0]
            ej.append(min(min(m[1] for m in h.dict()) for h in Dj) * (bs//(p-j)))
        Jy_exp = min(ej)
        cexp = n*factorial(p-1)
        f_ok = (Jy_exp == cexp)
        g_ok = ((n-p)*factorial(p-1) == cexp - bs) and (cexp % p != 0) and (((n-p)*factorial(p-1)) % p != 0)
        ok = all([a_ok,b_ok,c_ok,c1,c2,d_ok,e1,e2,e3,f_ok,g_ok]); allok &= ok
        print(f"p={p} n={n}: Sing=0 {a_ok}; d_x^(p)f=1 {b_ok}; chart y {c_ok}/{c1}; chart x no sing pts {c2}; "
              f"mu_mono={mu_mono} {d_ok}; arc w(f)=Np {e1}, mu_arc={mu_arc} {e2}, excludes d<p {e3}; "
              f"J#|x=0=(y^{cexp}) {f_ok}; exponents {g_ok}")
print("ALL CHECKS PASSED:", allok)

# (h) added 2026-10-24 for the general-p restructure: E_n is resolved by floor(n/p)
#     point blowups (always at the origin of the chart x = x'y), every blown-up point
#     carries an exponent > p, and the reduction exponent drops by p! each step:
#     c_{n-p} = c_n - p!.
allok2 = True
for p in [2,3,5,7]:
    for n in [m for m in range(p+1, 2*p+12) if m % p != 0][:6]:
        R.<x,y> = PolynomialRing(GF(p))
        g = x^p + y^n; steps = 0; exps = []
        while ordat(g, (0,0)) >= p:
            exps.append(min(sum(m) for m in (g - x^p).dict()))
            g = R(g(x*y, y) // y^p); steps += 1
        ok = (steps == n // p) and all(e > p for e in exps) and \
             all((e - p)*factorial(p-1) == e*factorial(p-1) - factorial(p) for e in exps)
        allok2 &= ok
        print(f"(h) p={p} n={n}: blowups={steps} (floor(n/p)={n//p}), exponents {exps}: {ok}")
print("(h) ALL CHECKS PASSED:", allok2)