# Sage script: computational checks for the adversarial review of
# Hironaka, "Resolution of Singularities in Positive Characteristics" (2017 preprint).
# Run with:  sage verify.sage   (tested with SageMath 10.8)
#
# Checks:
#  V1  E=(x^2+y^5,2) in char 2 (dim 2): singular loci along the resolution, tangent
#      cones, maximal-contact property of {x=0}, ambient reduction to {x=0},
#      and that Diff^(2)(f) contains a unit (so p^nega(E,-a) = O for all a>0).
#  V2  Narasimhan's example in char 2 (dim 4): the order-2 locus is the monomial
#      curve (t^32,t^7,t^19,t^15); its embedding dimension at 0 is 4; the first
#      derivatives of f have order >= 3 at 0 (so p(E,1) has no order-1 element at 0,
#      refuting the second sentence of Thm 6.1); at other points of the curve a
#      first derivative has order 1 (so the edge data jumps, Sing(E-check) = {0}).
#  V3  Hauser's kangaroo sequence (char 2, dim 3): x^2 + y^7 + y z^4, three point
#      blowups; orders stay 2; coefficient ideal w.r.t. {x=0} vs. the corrected
#      hypersurface x + y z^3 (1+y): residual order increases 2 -> 3.
#  V4  Edge data of the tangent cone (ridge) is non-increasing at near points for
#      random hypersurfaces of order 2 and 4 in char 2 (sanity check of Thm 7.1 at
#      the level of initial forms; NOT a check of the full p-algebra invariant).

import itertools, random
out = []
def say(*a):
    s = " ".join(str(x) for x in a)
    print(s); out.append(s)

def hasse(f, beta):
    """Hasse derivative d^(beta) f, beta a tuple of exponents, in f.parent()."""
    R = f.parent(); n = R.ngens(); 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(n):
            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 order_at(f, pt):
    """order of f at the point pt (coordinates in the base field or extension)."""
    R = f.parent()
    g = f(*[R.gen(i) + pt[i] for i in range(R.ngens())])
    if g == 0:
        return oo
    return min(sum(m) for m in g.dict().keys())

def initial_form(f, pt):
    R = f.parent()
    g = f(*[R.gen(i) + pt[i] for i in range(R.ngens())])
    d = min(sum(m) for m in g.dict().keys())
    return sum(c * R.monomial(*m) for m, c in g.dict().items() if sum(m) == d)

def singular_points_of_order(f, b, field):
    """closed points over `field` (finite) with ord_pt f >= b, brute force."""
    R = f.parent(); n = R.ngens()
    pts = []
    for pt in itertools.product(field, repeat=n):
        if order_at(f, pt) >= b:
            pts.append(pt)
    return pts

# ----------------------------------------------------------------------------------
say("="*78)
say("V1. E = (x^2 + y^5, 2), characteristic 2, ambient A^2")
F2 = GF(2)
R.<x,y> = PolynomialRing(F2)
f = x^2 + y^5
say("f =", f, "  ord_0 f =", order_at(f, (0,0)), "  in_0 f =", initial_form(f, (0,0)))
# singular locus = {ord >= 2} = {f = f_x = f_y = 0}
Jsing = R.ideal([f] + [hasse(f, b) for b in [(1,0),(0,1)]])
say("Sing(E) ideal (f, df):", Jsing.radical().gens(), " -> the origin only")
# Diff^(2) f contains a unit:
say("d_x^(2) f =", hasse(f, (2,0)), "  (a unit: so Diff^(2)(p(E,2)) = O, hence p^nega(E,-a)_0 = O for all a>0)")
# blow up the origin: chart y: x = x1*y, y = y
def blowup_point_chart(f, i):
    """chart where variable i is the exceptional parameter: x_j -> x_j * x_i (j != i)."""
    R = f.parent(); n = R.ngens(); g = R.gens()
    sub = [g[j] * g[i] if j != i else g[i] for j in range(n)]
    return f(*sub)
seq = [f]
cur = f
for step in range(3):
    # chart y (x = x' y): strict transform f' = f(x' y, y) / y^2
    tot = blowup_point_chart(cur, 1)
    cur = R(tot // y^2)
    assert tot == cur * y^2
    pts = singular_points_of_order(cur, 2, GF(4, 'a'))
    say(f"step {step+1}: chart y: f' = {cur};  ord_0 f' = {order_at(cur,(0,0))};  points of order>=2 over GF(4) on this chart: {pts}")
    # chart x (y = y' x): f(x, y' x)/x^2
    tot2 = blowup_point_chart(seq[-1], 0)
    other = R(tot2 // x^2)
    say(f"           chart x: f' = {other};  points of order>=2 over GF(4): {singular_points_of_order(other, 2, GF(4,'a'))}")
    seq.append(cur)
    if order_at(cur, (0,0)) < 2:
        say("   -> resolved (order < 2).")
        break
say("Observation: every point of order 2 of every transform lies on x'=0, the strict")
say("transform of {x=0}; the edge data (3-tuple (n, n-r, q)) of the p-algebra at the")
say("singular point is the same at every step (the only order-1 element of degree 1 is x),")
say("while the residual order ord(eps) drops 5 -> 3 -> 1.  The manuscript's Thm 16.6(2)")
say("asserts a strict lexicographic drop of the Inv-string at every point of pi^{-1}(nabla)\\nabla'.")
# ambient reduction to Y = {x = 0}: J# = sum_j (Diff^(j) J)^{b!/(b-j)}, b = 2, b! = 2
Diff1 = R.ideal([f, hasse(f,(1,0)), hasse(f,(0,1))])
Jsharp = R.ideal([f]) + Diff1^2
Ysub = R.ideal([x])
JY = (Jsharp + Ysub).quotient(R.ideal(1))  # just add x
JY = R.ideal(list(Jsharp.gens()) + [x])
say("J#|_{x=0} (generators mod x):", [g.reduce([x]) for g in Jsharp.gens() if g.reduce([x]) != 0], " -> F(Y) = (y^5, 2) on Y = A^1")

# ----------------------------------------------------------------------------------
say("="*78)
say("V2. Narasimhan's example, characteristic 2, ambient A^4")
S.<X,Y,Z,W> = PolynomialRing(F2)
fN = X^2 + Y*Z^3 + Z*W^3 + Y^7*W
say("f =", fN)
T.<t> = PolynomialRing(F2)
par = (t^32, t^7, t^19, t^15)
say("curve C: (x,y,z,w) = (t^32, t^7, t^19, t^15)")
say("f(C) =", fN(*par))
derivs = {v: hasse(fN, b) for v, b in zip("xyzw", [(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])}
for v, d in derivs.items():
    say(f"  d_{v} f = {d};   restricted to C: {d(*par)};   ord_0 = {order_at(d,(0,0,0,0))}")
crit = S.ideal([fN] + list(derivs.values()))
say("dimension of the critical scheme {f = df = 0}:", crit.dimension())
# ideal of C by elimination
U = PolynomialRing(F2, 'tt,X,Y,Z,W', order='lex')
tt, X_, Y_, Z_, W_ = U.gens()
IC = U.ideal([X_ - tt^32, Y_ - tt^7, Z_ - tt^19, W_ - tt^15]).elimination_ideal(tt)
IC_S = S.ideal([S(str(g)) for g in IC.gens()])
say("I(C) has", len(IC.gens()), "generators; minimal order of a generator at 0:",
    min(order_at(g,(0,0,0,0)) for g in IC_S.gens()), " (no order-1 element: C lies in no smooth hypersurface at 0)")
def is_min_gen(g, others):
    reach = [False]*(g+1); reach[0] = True
    for n in range(1, g+1):
        reach[n] = any(n-h >= 0 and reach[n-h] for h in others)
    return not reach[g]
sg = [32,7,19,15]
say("semigroup <32,7,19,15>: minimal generators:", [g for g in sg if is_min_gen(g,[h for h in sg if h != g])],
    " (all four => embedding dimension 4 at 0)")
# check I(C) = radical of critical ideal
say("radical(crit) == I(C):", crit.radical() == IC_S)
# ridge at a non-origin point of C: pick t0 in GF(2^k) with t0 != 0; use GF(4)
F4.<a> = GF(4)
S4 = S.change_ring(F4)
fN4 = S4(fN)
eta = tuple(F4(a)^e for e in (32, 7, 19, 15))
say("point eta on C over GF(4):", eta, " ord_eta f =", order_at(fN4, eta), " in_eta f =", initial_form(fN4, eta))
dz = S4(derivs['z'])
say("ord_eta(d_z f) =", order_at(dz, eta), " -> a degree-1 element of p(E,1) of order 1 at eta;",
    "so the edge data at eta is (4, <=2, ...) < (4,3,2) = edge data at 0.")

# ----------------------------------------------------------------------------------
say("="*78)
say("V3. Hauser's kangaroo sequence (arXiv:0811.4151, sect. G), characteristic 2, ambient A^3")
Q.<x,y,z> = PolynomialRing(F2)
f0 = x^2 + y^7 + y*z^4
say("f0 =", f0, "  Sing(f0,2) ideal radical:", Q.ideal([f0, hasse(f0,(0,1,0)), hasse(f0,(0,0,1)), hasse(f0,(1,0,0))]).radical().gens())
f1 = Q(f0(x*y, y, z*y) // y^2); say("f1 = f0(xy,y,zy)/y^2 =", f1)
f2 = Q(f1(x*z, y*z, z) // z^2); say("f2 = f1(xz,yz,z)/z^2 =", f2)
f3 = Q(f2(x*z, y*z + z, z) // z^2); say("f3 = f2(xz,yz+z,z)/z^2 =", f3)
for name, g in [("a0", f0), ("a1", f1), ("a2", f2), ("a3", f3)]:
    say(f"  ord at {name}: {order_at(g,(0,0,0))}, initial form: {initial_form(g,(0,0,0))}")
# residual orders w.r.t. the hypersurface x = 0 (coefficient ideal = eps = f - x^2 restricted), after
# removing exceptional monomials.
def residual_order_wrt_x(g, exc_vars):
    eps = g - x^2
    # divide out the exceptional monomial
    for v in exc_vars:
        while eps != 0 and all(m[Q.gens().index(v)] >= 1 for m in eps.dict().keys()):
            eps = Q(eps // v)
    return eps, order_at(eps, (0,0,0))
for name, g, exc in [("a0", f0, []), ("a1", f1, [y]), ("a2", f2, [y, z]), ("a3", f3, [z])]:
    eps, o = residual_order_wrt_x(g, exc)
    say(f"  {name}: coefficient part w.r.t. x=0 (exceptional monomial removed) = {eps}, residual order = {o}")
# at a3, the manuscript's LL-tail recipe: remove the square part of eps in the coordinates (y,z)
eps3 = f3 - x^2
sq = sum(c*Q.monomial(*m) for m, c in eps3.dict().items() if all(e % 2 == 0 for e in m))
nonsq = eps3 - sq
phi = sum(Q.monomial(*[e//2 for e in m]) for m, c in sq.dict().items())
say("  a3: square part of eps =", sq, " = (", phi, ")^2 ;  non-square part h =", nonsq)
say("  manuscript-style tail at a3: x + phi =", x + phi, " ; hypersurface {x + phi = 0}")
g3 = f3(x + phi, y, z)
say("  f3 in coordinates x -> x + phi:", g3)
eps_new, o_new = residual_order_wrt_x(g3, [z])
say("  residual order w.r.t. the new hypersurface:", o_new, " (was 2 at a2 w.r.t. {x=0}: INCREASE 2 -> 3)")
say("  the transform of the old tail {x=0} still contains Sing(f3,2) but has residual order",
    residual_order_wrt_x(f3,[z])[1], "which is NOT the residual order realised by the correct tail.")

# ----------------------------------------------------------------------------------
say("="*78)
say("V4. Ridge (edge data of the tangent cone) at near points, random hypersurfaces, char 2")
def ridge_data(F):
    """F homogeneous. Return (number of additive generators, sorted degrees) of the
    Diff-closure of (F) under Hasse derivatives of order < deg, computed via additive
    polynomials in the ideal."""
    R = F.parent(); n = R.ngens(); K = R.base_ring(); p = K.characteristic()
    gens = [F]; I = R.ideal(gens)
    changed = True
    while changed:
        changed = False
        new = []
        for g in list(gens):
            d = g.degree()
            for beta in itertools.product(range(d), repeat=n):
                if 0 < sum(beta) < d:
                    h = hasse(g, beta)
                    if h != 0 and h not in I:
                        new.append(h)
        if new:
            gens = gens + new
            I = R.ideal(gens); changed = True
    # additive polynomials in I by degree p^k
    dims = []
    k = 0
    D = F.degree()
    while p^k <= D:
        q = p^k
        # linear algebra: which combos sum a_i x_i^q lie in I
        GB = I.groebner_basis()
        vecs = []
        monos = set()
        reds = []
        for i in range(n):
            r = (R.gen(i)^q).reduce(GB)
            reds.append(r)
            monos |= set(r.dict().keys())
        monos = sorted(monos)
        M = matrix(K, n, len(monos), lambda i, j: reds[i].dict().get(monos[j], 0))
        dims.append(M.left_kernel().dimension())
        k += 1
    # multiplicities: m_k = dims[k] - dims[k-1]
    degs = []
    prev = 0
    for k, dk in enumerate(dims):
        for _ in range(dk - prev):
            degs.append(p^k)
        prev = dk
    return (len(degs), degs)

random.seed(int(1))
R3.<x,y,z> = PolynomialRing(F2)
viol = 0; tested = 0
for trial in range(60):
    b = random.choice([2, 4])
    # random f with ord_0 = b, initial form a p-th power-ish: force in_0 = x^b plus random higher terms
    f = x^b
    for _ in range(random.randint(1, 5)):
        e = [random.randint(0, 3) for _ in range(3)]
        if sum(e) > b:
            f += R3.monomial(*e)
    if order_at(f, (0,0,0)) != b:
        continue
    inv0 = ridge_data(initial_form(f, (0,0,0)))
    # blow up the origin, all three charts, points over GF(4) on exceptional divisor with order b
    F4r = GF(4, 'a'); R4 = R3.change_ring(F4r); f4 = R4(f)
    for i in range(3):
        g = R4.gens(); sub = [g[j]*g[i] if j != i else g[i] for j in range(3)]
        tot = f4(*sub); fp = R4(tot // g[i]^b)
        for pt in itertools.product(F4r, repeat=3):
            if pt[i] != 0:
                continue
            if order_at(fp, pt) == b:
                inv1 = ridge_data(initial_form(fp, pt))
                tested += 1
                # compare (n - r, degs) lexicographically: inv1 should be <= inv0
                key0 = (3 - inv0[0], tuple(inv0[1])); key1 = (3 - inv1[0], tuple(inv1[1]))
                if key1 > key0:
                    viol += 1
                    say("  VIOLATION?", f, "->", fp, "at", pt, inv0, inv1)
say(f"  near points tested: {tested}; lexicographic increases of (n-r, q) found: {viol}")

open("verify-output.txt", "w").write("\n".join(out) + "\n")
