"""Exact finite Descartes-generator orbits; no floating point or external code.

Bourgain–Kontorovich arXiv:1205.4416v2, Section 2.1: Gamma is generated
by S1*S2, S2*S3, S3*S4. Orbit-stabilizer identifies the vector orbit with
Gamma / Stab(v0 mod p). Finite data alone do not establish an asymptotic.
"""
from collections import Counter,deque
from itertools import product
from math import prod
import json

def reflect(v,i,p=None):
    x=list(v);x[i]=2*sum(v)-3*v[i]
    return tuple(a%p for a in x) if p else tuple(x)

def act(v,pair,p=None):
    return reflect(reflect(v,pair[1],p),pair[0],p)

def form(v):return 2*sum(x*x for x in v)-sum(v)**2

def orbit(v,p,pairs):
    v=tuple(x%p for x in v);seen={v};todo=deque([v])
    while todo:
        u=todo.popleft()
        for pair in pairs:
            w=act(u,pair,p)
            if w not in seen:seen.add(w);todo.append(w)
    return seen

pairs=[(0,1),(1,2),(2,3)]
v0=(-11,21,24,28)
assert form(v0)==0
rows=[]
for p in [5,7,11,13,17,19]:
    orb=orbit(v0,p,pairs)
    null={v for v in product(range(p),repeat=4) if any(v) and form(v)%p==0}
    assert orb==null,(p,len(orb),len(null))
    fibers=Counter(v[0] for v in orb)
    legendre=1 if p%4==1 else -1
    formula=(p-1)*(p*p+(legendre+1)*p+1)
    assert len(orb)==formula
    rows.append({'p':p,'vector_stabilizer_index':len(orb),'nonzero_null_vectors':len(null),
        'projective_orbit_size':len(orb)//(p-1),'p_squared':p*p,
        'first_coordinate_fiber_counts':dict(sorted(fibers.items())),
        'finite_field_formula':formula,'comparison_only':'Finite checks; the all-prime claim requires a separate proof.'})
# The full SO contains -I, but all generators act on the same finite group;
# group closure at p=7 directly distinguishes full SO from its spin image.
p=7
I=tuple(1 if i==j else 0 for i in range(4) for j in range(4))
minusI=tuple((-x)%p for x in I)
def matmul(a,b):return tuple(sum(a[4*i+k]*b[4*k+j] for k in range(4))%p for i in range(4) for j in range(4))
matrices=[]
for pair in pairs:
    cols=[act(tuple(1 if i==j else 0 for i in range(4)),pair,p) for j in range(4)]
    matrices.append(tuple(cols[j][i] for i in range(4) for j in range(4)))
group={I};todo=deque([I])
while todo:
    a=todo.popleft()
    for b in matrices:
        c=matmul(a,b)
        if c not in group:group.add(c);todo.append(c)
assert minusI not in group
# The exact closure, rather than a guessed orthogonal-group order, determines
# the image cardinality. Adjoining central -I gives twice as many isometries.
assert all(matmul(minusI,a) not in group for a in group)
# -I preserves F and has determinant +1 in dimension four.
# Thus it lies in SO_F(F7), independently of any classification theorem.
result={'arithmetic':'exact Python integers and modular residues','root':v0,'generators':'S1S2, S2S3, S3S4 (rightmost reflection acts first)',
    'orbit_checks':rows,'group_check_mod_7':{'image_order':len(group),'minus_identity_in_image':minusI in group,
    'SO_order_lower_bound_by_adjoining_central_minus_I':2*len(group),'minus_identity_preserves_form':True,'minus_identity_determinant':1,
    'conclusion':'Image has 58800 elements; adjoining the missing central -I yields 117600 distinct form-preserving determinant-one matrices. Thus full SO_F(F7) has strictly larger order than the printed quotient, independently of a classification theorem.'}}
print(json.dumps(result,indent=2))