# Check that Int(sigma_{w0}) o C is the identity on so(8) (split D4, J antidiagonal), C(X) = -X^T.
# Tits lift sigma_a = x_a(1) x_{-a}(-1) x_a(1); in the vector representation x_a(c) = exp(c e_a) = 1 + c e_a.
D = 8
def zero(): return [[0]*D for _ in range(D)]
def E(i, j):
    m = zero(); m[i-1][j-1] = 1; return m
def add(a, b, s=1): return [[a[i][j] + s*b[i][j] for j in range(D)] for i in range(D)]
def mul(a, b): return [[sum(a[i][k]*b[k][j] for k in range(D)) for j in range(D)] for i in range(D)]
def tr(a): return [[a[j][i] for j in range(D)] for i in range(D)]
I = [[int(i == j) for j in range(D)] for i in range(D)]
J = [[int(i + j == D - 1) for j in range(D)] for i in range(D)]
dual = lambda i: D + 1 - i
# simple roots of D4: e1-e2, e2-e3, e3-e4, e3+e4 (Bourbaki D4, centre = node 2)
e = [add(E(1,2), E(dual(2),dual(1)), -1), add(E(2,3), E(dual(3),dual(2)), -1),
     add(E(3,4), E(dual(4),dual(3)), -1), add(E(3,dual(4)), E(4,dual(3)), -1)]
f = [tr(x) for x in e]
for x in e + f:   # in so(8)_J and nilpotent of square 0
    assert add(mul(tr(x), J), mul(J, x)) == zero() and mul(x, x) == zero()
C = lambda X: [[-v for v in row] for row in tr(X)]
sig = [mul(mul(add(I, e[k]), add(I, f[k], -1)), add(I, e[k])) for k in range(4)]
# reduced word of w0 for D4 via the root action on coefficient vectors
Ad4 = [[2,-1,0,0],[-1,2,-1,-1],[0,-1,2,0],[0,-1,0,2]]
def refl(b, i):
    c = sum(b[j]*Ad4[j][i] for j in range(4)); return tuple(b[j] - (c if j == i else 0) for j in range(4))
simple = [tuple(int(j == i) for j in range(4)) for i in range(4)]
pos = set(simple); fr = list(simple)
while fr:
    new = []
    for b in fr:
        for i in range(4):
            r = refl(b, i)
            if all(t >= 0 for t in r) and r not in pos: pos.add(r); new.append(r)
    fr = new
assert len(pos) == 12
def act(word, b):   # apply s_{word[0]} ... s_{word[-1]} to b (rightmost first)
    for i in reversed(word): b = refl(b, i)
    return b
word = []
while True:
    js = [j for j in range(4) if all(t >= 0 for t in act(word, simple[j]))]
    if not js: break
    word.append(js[0])
assert len(word) == 12
s0 = I
for k in word: s0 = mul(s0, sig[k])
s0inv = tr(s0)   # orthogonal? check
assert mul(s0, s0inv) == I or True
# s0 lies in O(8)_J: s0^T J s0 = J ; inverse = J s0^T J
assert mul(mul(tr(s0), J), s0) == J
s0inv = mul(mul(J, tr(s0)), J)
assert mul(s0, s0inv) == I
ok = all(mul(mul(s0, C(x)), s0inv) == x for x in e + f)
print("reduced word of w0(D4), 1-based nodes:", [k+1 for k in word])
print("Ad(sigma_w0) o C fixes all e_i, f_i:", ok)
print("sigma_w0^2 = 1 in SO(8):", mul(s0, s0) == I)