#!/usr/bin/python3

# see also: https://oeis.org/A093199

from ortools.constraint_solver import pywrapcp

def main():
  solver = pywrapcp.Solver("...")

  cells=[[solver.IntVar(0, 9, "%d_%d" % (r, c)) for r in range(4)] for c in range(4)]
  sum_=solver.IntVar(0, 9*4, "sum_")

  #print (cells)

  for r in range(4):
      solver.Add(cells[r][0]+cells[r][1]+cells[r][2]+cells[r][3]==sum_)
  for c in range(4):
      solver.Add(cells[0][c]+cells[1][c]+cells[2][c]+cells[3][c]==sum_)
  solver.Add(cells[0][0]+cells[1][1]+cells[2][2]+cells[3][3]==sum_)
  solver.Add(cells[3][0]+cells[2][1]+cells[1][2]+cells[0][3]==sum_)

  solution = solver.Assignment()

  # flatten
  x=[]
  for r in range(4):
      for c in range(4):
          x.append(cells[r][c])

  db = solver.Phase(x, solver.CHOOSE_MIN_SIZE_LOWEST_MAX,
                    solver.ASSIGN_MIN_VALUE)

  solver.NewSearch(db)
  solutions = 0
  while solver.NextSolution():
      #for r in range(4):
      #    for c in range(4):
      #        print ("%d " % cells[r][c].Value(), end="")
      #    print ("")
      #print ("")
      solutions+=1
  solver.EndSearch()
  return solutions

solutions = main()
print (f"{solutions=}")