
PROG

(SageMath)
def msos_search(F, single=False):
squares = {x^2 for x in F}
MSOS = []
E = 0
for A, I in Subsets(squares, 2):
if A + I != 2*E: continue
C, G = 1, 1
B = 3*E  A  C
D = 3*E  A  G
F = 3*E  C  I
H = 3*E  G  I
if len(squares & {B, D, F, H}) < 4: continue
if len({A, B, C, D, E, F, G, H, I}) < 9: continue
if single: return [A, B, C, D, E, F, G, H, I]
MSOS.append([A, B, C, D, E, F, G, H, I])
E = 1
sequences = []
for A, I in Subsets(squares, 2):
if A + I != 2*E: continue
for C, G in sequences:
B = 3*E  A  C
D = 3*E  A  G
F = 3*E  C  I
H = 3*E  G  I
if len(squares & {B, D, F, H}) < 4: continue
if len({A, B, C, D, E, F, G, H, I}) < 9: continue
if single: return [A, B, C, D, E, F, G, H, I]
MSOS.append([A, B, C, D, E, F, G, H, I])
sequences.append((A, I))
return MSOS
for q in range(3, 500):
if q % 2 == 0: continue
if len(factor(q)) > 1: continue
print q, msos_search(GF(q), single=True)
