"""Check: L(S) = (1/18) * sum_W a_W b_W c_W - C(n,3)/2 for S subset of F_3^4,
where W ranges over the 40 parallel classes of hyperplanes and (a,b,c) are the
slice sizes. Also the 3D analogue L = (1/6) sum_P abc - C(n,3)/2 over 13 plane classes."""
import itertools, random
from math import comb

def points(d):
    return list(itertools.product(range(3), repeat=d))

def lines_count(S):
    Sset = set(S)
    cnt = 0
    for x, y in itertools.combinations(S, 2):
        z = tuple((-a - b) % 3 for a, b in zip(x, y))
        if z in Sset:
            cnt += 1
    assert cnt % 3 == 0
    return cnt // 3

def classes(d):
    # nonzero xi up to sign
    seen = set(); out = []
    for xi in points(d):
        if any(xi):
            m = tuple((-a) % 3 for a in xi)
            if m in seen: continue
            seen.add(xi); out.append(xi)
    return out

def profile(S, xi):
    k = [0, 0, 0]
    for x in S:
        k[sum(a * b for a, b in zip(x, xi)) % 3] += 1
    return k

def check(d, trials):
    P = points(d); W = classes(d)
    assert len(W) == (3**d - 1) // 2
    denom = {4: 18, 3: 6}[d]
    for _ in range(trials):
        n = random.randint(0, 3**d)
        S = random.sample(P, n)
        L = lines_count(S)
        sabc = sum(eval('*'.join(map(str, profile(S, xi)))) for xi in W)
        rhs = sabc / denom - comb(n, 3) / 2
        assert abs(L - rhs) < 1e-9, (d, n, L, rhs)
        # pair identity
        se2 = 0
        for xi in W:
            a, b, c = profile(S, xi); se2 += a*b + b*c + c*a
        assert se2 == {4: 27, 3: 9}[d] * comb(n, 2)
    print(f"d={d}: identity holds on {trials} random sets")

random.seed(1)
check(4, 300)
check(3, 300)