How to rearrange numbers in a 4x4 grid to all sum up to 0

so I need to code a solution to see if by rearranging the numbers of a 4x4 grid to see if the columns, rows and diaganals all equal 0. I've been trying to solve this question for a bit but its not working.

I looked online and the only answer I could see was the brute force method which would make every possible method, the only problem is that there are 16! combinations aka 20 trillion combinations so I tried to flatten it and that still didn't work.

Have a look at the code please:

import pandas as pd
import numpy as np
from itertools import permutations, islice

# Initial grid
grid = pd.DataFrame({
    'A': [5,-5,3,-2],
    'B': [1,-1,9,-10],
    'C': [-7,-4,7,4],
    'D': [-3,-8,2,8]
})

# there are like 16! possible combinations (20 trillion)
def is_valid_solution(df):
    # Check rows and columns sum to 0
    if not (df.sum(axis=0) == 0).all() or not (df.sum(axis=1) == 0).all():
        return False
    # Check diagonals sum to 0
    if df.values.trace() != 0 or np.fliplr(df.values).trace() != 0:
        return False
    return True

def rearrange_grid_bruteforce(df):
    values = df.values.flatten()
    # Use islice to limit the number of permutations evaluated
    for perm in islice(permutations(values), 100000):  # Just as an example, not feasible for the entire set
        new_grid = np.array(perm).reshape(df.shape)
        new_df = pd.DataFrame(new_grid, columns=df.columns)
        if is_valid_solution(new_df):
            return new_df
    return None

# Attempt to find a solution
# Note: This is for demonstration and will not compute all permutations due to computational constraints
solution = rearrange_grid_bruteforce(grid)

print('Working on it....')

if solution is not None:
    print("A solution was found:")
    print(solution)
else:
    print("No solution found within the limited permutation set.")

Solution

You can solve this problem with binary integer linear programming. I adapted this test of SciPy's milp function to solve the case in which the numbers are given and the sum M=0. As indicated in the comments, your particular case is infeasible, but if the first number of the second row is changed to 2, a solution can be found in a few seconds.

import numpy as np
from scipy.optimize import milp

np.random.seed(0)

# Original code
# M = n * (n ** 2 + 1) / 2
# numbers = np.arange(n ** 4) // n ** 2 + 1
# numbers = numbers.reshape(n ** 2, n, n)

# Adapted to your problem
M = 0
numbers = [5, -5, 3, -2, 1, -1, 9, -10, -7, -4, 7, 4, -3, -8, 2, 8]
numbers[4] = 2  # change to make feasible
n = int(np.sqrt(len(numbers)))
numbers = np.reshape(numbers, (n**2, 1, 1))
numbers = np.tile(numbers, (1, n, n))

zeros = np.zeros((n ** 2, n, n))

A_list = []
b_list = []

# Rule 1: use every number exactly once
for i in range(n ** 2):
    A_row = zeros.copy()
    A_row[i, :, :] = 1
    A_list.append(A_row.flatten())
    b_list.append(1)

# Rule 2: Only one number per square
for i in range(n):
    for j in range(n):
        A_row = zeros.copy()
        A_row[:, i, j] = 1
        A_list.append(A_row.flatten())
        b_list.append(1)

# Rule 3: sum of rows is M
for i in range(n):
    A_row = zeros.copy()
    A_row[:, i, :] = numbers[:, i, :]
    A_list.append(A_row.flatten())
    b_list.append(M)

# Rule 4: sum of columns is M
for i in range(n):
    A_row = zeros.copy()
    A_row[:, :, i] = numbers[:, :, i]
    A_list.append(A_row.flatten())
    b_list.append(M)

# Rule 5: sum of diagonals is M
A_row = zeros.copy()
A_row[:, range(n), range(n)] = numbers[:, range(n), range(n)]
A_list.append(A_row.flatten())
b_list.append(M)
A_row = zeros.copy()
A_row[:, range(n), range(-1, -n - 1, -1)] = \
    numbers[:, range(n), range(-1, -n - 1, -1)]
A_list.append(A_row.flatten())
b_list.append(M)

A = np.array(np.vstack(A_list), dtype=float)
b = np.array(b_list, dtype=float)
c = np.random.rand(A.shape[1])  # arbitrary; we just need to satisfy constraints

res = milp(c=c * 0, constraints=(A, b, b), bounds=(0, 1), integrality=1)

x = np.round(res.x)
s = (numbers.flatten() * x).reshape(n ** 2, n, n)
square = np.sum(s, axis=0)
np.testing.assert_allclose(square.sum(axis=0), M)
np.testing.assert_allclose(square.sum(axis=1), M)
np.testing.assert_allclose(np.diag(square).sum(), M)
np.testing.assert_allclose(np.diag(square[:, ::-1]).sum(), M)
np.testing.assert_allclose(np.sort(square.ravel()), np.sort(square.ravel()))

square
# array([[ -3.,   2.,   9.,  -8.],
#        [  2.,  -5.,  -4.,   7.],
#        [ -7.,   4.,   5.,  -2.],
#        [  8.,  -1., -10.,   3.]])