main_ortools.py

from ortools.sat.python import cp_model

N = 12

model = cp_model.CpModel()

f = [[model.NewIntVar(0, N - 1, f'f{i}{j}') for j in range(N)] for i in range(N)]
g = [[model.NewIntVar(0, N - 1, f'g{i}{j}') for j in range(N)] for i in range(N)]
s = [[model.NewIntVar(0, N * N - 1, f's{i}{j}') for j in range(N)] for i in range(N)]

for i in range(N):
    model.Add(f[0][i] == i)
    model.Add(f[i][0] == i)
    model.Add(g[0][i] == i)

for i in range(N):
    for j in range(N):
        model.Add(s[i][j] == f[i][j] + N * g[i][j])

for i in range(N):
    for x in [f, g]:
        model.AddAllDifferent([x[i][j] for j in range(N)])
        model.AddAllDifferent([x[j][i] for j in range(N)])

model.AddAllDifferent([s[i][j] for i in range(4) for j in range(4)])

solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL:
    print("OPTIMAL")
elif status == cp_model.FEASIBLE:
    print("FEASIBLE")
elif status == cp_model.INFEASIBLE:
    print("INFEASIBLE")
    exit()
else:
    print(status)

print([[solver.Value(x) for x in line] for line in f])
print([[solver.Value(x) for x in line] for line in g])
print([[solver.Value(x) for x in line] for line in s])

for i in range(N):
    print(' '.join(f'{solver.Value(f[i][j])}{solver.Value(g[i][j])}' for j in range(N)))