"""Check the rank of the 144 cyclic Freudenthal quadrics over Q.

Run with Python 3; no third-party packages are required. The eight
relations are transcribed from the displayed defining equations.
"""
from fractions import Fraction
from pathlib import Path
import json

VARIABLES = [f'{kind}{i}' for kind in ('a', 'x', 'y') for i in range(18)] + ['z1', 'z2']
INDEX = {v: i for i, v in enumerate(VARIABLES)}
# Each entry is (coefficient, first variable, second variable).
RELATIONS = [
 [(1,'x1','x2'),(-1,'a3','x17'),(-1,'a2','y0'),(-1,'a1','y3'),(-1,'a0','x4'),(-1,'a17','a4')],
 [(1,'y14','y1'),(-1,'x15','y4'),(1,'a0','x10'),(1,'a15','x5'),(-1,'x0','y11'),(-1,'x13','x2')],
 [(1,'y2','y1'),(-1,'a2','y12'),(-1,'y5','x0'),(-1,'y16','x3'),(-1,'a1','y9'),(-1,'a5','a16')],
 [(1,'x4','y1'),(-1,'x2','y5'),(-1,'a5','x17'),(-1,'a4','y16'),(1,'a2','y8'),(1,'a1','x7')],
 [(1,'x2','x6'),(-1,'x4','y4'),(1,'a4','z2'),(1,'a5','y0'),(1,'a3','y8'),(-1,'a1','a7')],
 [(1,'y1','y3'),(-1,'x17','x5'),(-1,'x2','z1'),(-1,'a0','x7'),(-1,'a4','x15'),(1,'a2','y11')],
 [(1,'y4','y16'),(-1,'y1','z2'),(1,'y8','x0'),(1,'y12','x2'),(1,'a1','x10'),(-1,'a5','a15')],
 [(1,'z1','z2'),(-1,'y3','y12'),(-1,'y7','y16'),(-1,'y8','y17'),(-1,'x1','x10'),(-1,'x3','x12'),(-1,'x5','x14'),(-1,'a0','a9'),(-1,'a2','a11'),(1,'a3','a12'),(-1,'a4','a13'),(-1,'a6','a15')]
]

def shift(v, k):
    if v[0] == 'z':
        return 'z'+str(3-int(v[1])) if k % 2 else v
    return v[0]+str((int(v[1:])+k) % 18)

def rows():
    out = []
    for relation in RELATIONS:
        for k in range(18):
            row = {}
            for coefficient, u, v in relation:
                pair = tuple(sorted((INDEX[shift(u,k)], INDEX[shift(v,k)])))
                row[pair] = row.get(pair, Fraction(0)) + coefficient
            out.append({c:a for c,a in row.items() if a})
    return out

def exact_rank(matrix):
    echelon = {}
    for original in matrix:
        row = original.copy()
        while row:
            pivot = min(row)
            if pivot not in echelon:
                coefficient = row[pivot]
                echelon[pivot] = {c:a/coefficient for c,a in row.items()}
                break
            coefficient = row[pivot]
            for c,a in echelon[pivot].items():
                row[c] = row.get(c, Fraction(0))-coefficient*a
                if not row[c]:
                    del row[c]
    return len(echelon)

def main():
    matrix = rows()
    rank = exact_rank(matrix)
    assert len(VARIABLES)==56
    assert len(matrix)==144
    assert rank==133, f'Unexpected rank: {rank}'
    report = {'variables':56,'cyclic_quadrics':144,'rank_over_Q':rank,
              'nullity':144-rank,'method':'Exact sparse Gaussian elimination over Fraction',
              'scope':'Linear rank of the printed quadrics, not source-ring generation.'}
    Path(__file__).with_name('freudenthal_quadrics_report.json').write_text(json.dumps(report,indent=2)+'\n')
    print(json.dumps(report,indent=2))

if __name__ == '__main__':
    main()