import sys, os; sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
# (c) theta0 = Int(sigma_L) o C on E6 restricted to the black D4: is Ad(sigma_{w0}) o C the identity on so(8)?
# Pinned split so(8): basis e1..e4,f1..f4 (indices 0..7), form J = [[0,I],[I,0]].
N = 8
def zero(): return [[0]*N for _ in range(N)]
def E(i, j): m = zero(); m[i][j] = 1; return m
def add(a, b, s=1): return [[a[i][j] + s*b[i][j] for j in range(N)] for i in range(N)]
def mul(a, b): return [[sum(a[i][k]*b[k][j] for k in range(N)) for j in range(N)] for i in range(N)]
def tr(a): return [[a[j][i] for j in range(N)] for i in range(N)]
I8 = [[int(i == j) for j in range(N)] for i in range(N)]
J = zero()
for i in range(4): J[i][4+i] = 1; J[4+i][i] = 1
def in_so(X): return add(mul(tr(X), J), mul(J, X)) == zero()
# simple roots of D4: a1=e1-e2, a2=e2-e3 (centre), a3=e3-e4, a4=e3+e4
def Eroot_minus(i, j): return add(E(i, j), E(4+j, 4+i), -1)        # eps_i - eps_j
def Eroot_plus(i, j): return add(E(i, 4+j), E(j, 4+i), -1)          # eps_i + eps_j
e = [Eroot_minus(0,1), Eroot_minus(1,2), Eroot_minus(2,3), Eroot_plus(2,3)]
f = [tr(x) for x in e]
assert all(in_so(x) for x in e + f)
h = [add(mul(e[i], f[i]), mul(f[i], e[i]), -1) for i in range(4)]
# Cartan matrix check: [h_i, e_j] = A_ij e_j
def br(a, b): return add(mul(a, b), mul(b, a), -1)
Acart = [[None]*4 for _ in range(4)]
for i in range(4):
    for j in range(4):
        c = br(h[i], e[j])
        for k in (2, 0, -1, -2, 1):
            if c == [[k*x for x in row] for row in e[j]]: Acart[i][j] = k
print('Cartan matrix <alpha_j, h_i>:', Acart)
def expnil(x):   # x^2 = 0 for these root vectors
    assert mul(x, x) == zero(); return add(I8, x)
sig = [mul(mul(expnil(e[i]), expnil([[-v for v in r] for r in f[i]])), expnil(e[i])) for i in range(4)]
# reduced word of w0 via root action: simple roots in basis eps
roots_s = [(1,-1,0,0), (0,1,-1,0), (0,0,1,-1), (0,0,1,1)]
def refl(i, v):
    a = roots_s[i]; c = sum(x*y for x, y in zip(a, v)); return tuple(x - c*y for x, y in zip(v, a))
def wact(word, v):
    for i in reversed(word): v = refl(i, v)
    return v
def positive(v):
    for x in v:
        if x != 0: return x > 0
    return True
word = []
while True:
    ext = [i for i in range(4) if positive(wact(word, roots_s[i]))]
    if not ext: break
    word.append(ext[0])
print('reduced word of w0 (0-based):', word, 'length', len(word))
assert all(wact(word, r) == tuple(-x for x in r) for r in roots_s)   # w0 = -1
S = I8
for i in word: S = mul(S, sig[i])
Sinv = tr(S)            # S orthogonal? check
assert mul(S, Sinv) == I8 or True
def inv(M):
    # integer inverse via J-orthogonality: M^{-1} = J M^T J for M in O(J)
    return mul(mul(J, tr(M)), J)
assert mul(S, inv(S)) == I8
print('sigma_{w0}^2 == I in SO(8):', mul(S, S) == I8)
def C(X): return [[-v for v in r] for r in tr(X)]
ok = all(mul(mul(S, C(e[i])), inv(S)) == e[i] for i in range(4)) and all(mul(mul(S, C(f[i])), inv(S)) == f[i] for i in range(4))
print('Ad(sigma_{w0}) o C fixes all e_i, f_i (hence is the identity on so(8)):', ok)