residues = sorted(set((i*i) % 31 for i in range(1, 31)))
nonresidues = sorted(set(range(1, 31)).difference(residues))

omega = 2
def mul(a,h,b,k):
    if a == b:
        return (a, (4*h+k) % 7)
    if (a - b) % 31 in residues:
        return ((3*a-2*b) % 31, k % 7)
    return ((3*a-2*b) % 31, (h - omega * (k - h)) % 7)

magma217 = [[7 * mul(a, h, b, k)[0] + mul(a, h, b, k)[1]
             for b in range(31) for k in range(7)]
            for a in range(31) for h in range(7)]
#print(all(i == magma217[j][magma217[i][magma217[magma217[j][i]][j]]] for i in range(217) for j in range(217)))
#print(all(len(set(magma217[i])) == 217 for i in range(217)))
#print(all(len(set(magma217[j][i] for j in range(217))) == 127 for i in range(217)))

print(magma217)