How to find a Magic Square of Squares (n^2)?

Question

I am trying to find a magic square whose elements are distinct squares of natural numbers.

I have 200 solutions of 4 numbers, for each row or column of the magic square. But how would I combine them to 16 numbers whose sum is 8515, one of possible sums in a magic square of squares?

And how would we make a generalized program that works for any given numbers, and n x n.

import math
import time
from collections import Counter

def find_solutions(N):
    solutions = set()
    sqrt_N = int(math.isqrt(N))  # Calculate the square root of N
    for a in range(sqrt_N + 1):
        a_squared = a**2
        for b in range(sqrt_N + 1):  # Start from 0
            b_squared = b**2
            for c in range(sqrt_N + 1):  # Start from 0
                c_squared = c**2
                remaining = N - a_squared - b_squared - c_squared

                if remaining < 0:
                    break  # Stop if remaining is negative

                d = int(math.isqrt(remaining))
                d_squared = d**2

                if a_squared + b_squared + c_squared + d_squared == N and len({a, b, c, d}) == 4:
                    solutions.add(tuple(sorted([a, b, c, d])))

    return solutions

# Specify your value of N
N = 8515  # Change this to your desired value

# Print the result and execution time
solution_count = len(solutions)

if solution_count == 0:
    print(f"No solutions found for a^2 + b^2 + c^2 + d^2 = {N}")
else:
    print(f"Number of unique solutions for a^2 + b^2 + c^2 + d^2 = {N}: {solution_count}")

    # Calculate total occurrences for each number in the list of solutions
    total_counts = Counter([num for sublist in solutions for num in sublist])

    # Calculate percentages for each number
    total_solutions = sum(total_counts.values())
    number_percentages = {num: (count / total_solutions) * 100 for num, count in total_counts.items()}

    # Sort solutions by total percentage in descending order
    sorted_solutions = sorted(solutions, key=lambda solution: sum(number_percentages[num] for num in solution), reverse=True)

    # Print the results
    print("\nSolutions in Descending Order of Total Percentage:")
    print("{:<30} {:<10}".format("[a, b, c, d]", "Total Percentage"))
    for solution in sorted_solutions:
        total_percentage = sum(number_percentages[num] for num in solution)
        print("{:<30} {:<10.5f}%".format(str(solution), total_percentage))

Well i have been trying brute force, i have been trying different search algortihms but nothing works for me, and i just need to be over it. I have been trying graph theory, but it takes too much time. I even went as far as to imagine each values in the magic square of square as function cause it is like qudratic function f(x) = a^2 + k . And then somehow optimize it. But, it didn't work for me. I am requesting help from the community! Thanks.

UPDATE_______

Thanks for all you guys questions. Now, i have updated it slightly, and now i want to find a solutions which is "faster."

Here was my idea, but i am stuck on as to how to even match it so it solves the magic square of square :

import math
import time
from collections import Counter

def find_solutions(N):
    solutions = set()
    sqrt_N = int(math.isqrt(N))  # Calculate the square root of N
    for a in range(sqrt_N + 1):
        a_squared = a**2
        for b in range(sqrt_N + 1):  # Start from 0
            b_squared = b**2
            for c in range(sqrt_N + 1):  # Start from 0
                c_squared = c**2
                remaining = N - a_squared - b_squared - c_squared

                if remaining < 0:
                    break  # Stop if remaining is negative

                d = int(math.isqrt(remaining))
                d_squared = d**2

                if a_squared + b_squared + c_squared + d_squared == N and len({a, b, c, d}) == 4:
                    solutions.add(tuple(sorted([a, b, c, d])))

    return solutions

# Specify your value of N
N = 8515  # Change this to your desired value

# Print the result and execution time
solution_count = len(solutions)

if solution_count == 0:
    print(f"No solutions found for a^2 + b^2 + c^2 + d^2 = {N}")
else:
    print(f"Number of unique solutions for a^2 + b^2 + c^2 + d^2 = {N}: {solution_count}")

    # Calculate total occurrences for each number in the list of solutions
    total_counts = Counter([num for sublist in solutions for num in sublist])

    # Calculate percentages for each number
    total_solutions = sum(total_counts.values())
    number_percentages = {num: (count / total_solutions) * 100 for num, count in total_counts.items()}

    # Sort solutions by total percentage in descending order
    sorted_solutions = sorted(solutions, key=lambda solution: sum(number_percentages[num] for num in solution), reverse=True)

    # Print the results
    print("\nSolutions in Descending Order of Total Percentage:")
    print("{:<30} {:<10}".format("[a, b, c, d]", "Total Percentage"))
    for solution in sorted_solutions:
        total_percentage = sum(number_percentages[num] for num in solution)
        print("{:<30} {:<10.5f}%".format(str(solution), total_percentage))

Answer 1

First of all, when finding solutions for 4 numbers (a, b, c, d), you can introduce the constraint a < b < c < d into the search by adjusting the ranges. You don't need to check if some of the numbers are equal; you don't need to sort them, and you can arrange the solutions as a list and not a set (list is faster to append to).

def find_solutions(N):
    solutions = []
    sqrt_N = int(math.isqrt(N))  # Calculate the square root of N
    for a in range(0, sqrt_N):
        a_squared = a**2
        for b in range(a+1, sqrt_N):  # Start from a+1, not from 0
            b_squared = b**2
            for c in range(b+1, sqrt_N):  # Start from b+1, not from 0
                c_squared = c**2
                remaining = N - a_squared - b_squared - c_squared

                if remaining <= c_squared:
                    break  # Stop if remaining is too small

                d = int(math.isqrt(remaining))
                d_squared = d**2

                if a_squared + b_squared + c_squared + d_squared == N:
                    solutions.append((a, b, c, d))

    return solutions

To find a magic square, you should combine 4 row solutions, and check that the columns don't violate the magic-square conditions. A brute-force search could be something like this:

for a, b, c, d in solutions:
    for e, f, g, h in solutions:
        for i, j, k, l in solutions:
            for m, n, o, p in solutions:
                if (whatever):
                    # print the solution

To improve the performance of the search, you can do the iteration in a different order. There are a lot of constraints on the numbers in magic squares; the guideline is to check the constraints as close to the outer loop as possible. Changing the order of iteration makes it possible to apply constraints early.

I noticed that for each pair of adjacent numbers in a row of a magic square, the choice for the rest of them is very limited. For example, if there are 80 and 25 in a column or row, the other numbers are either 11 and 37 or 23 and 31 — only 4 possibilities (taking order into account)!

To take advantage of this, make a dictionary which says which possibilities exist for each pair of numbers.

poss = {}
for s in solutions:
    for a, b, c, d in itertools.permutations(s):
        if (a, b) not in poss:
            poss[(a, b)] = []
        poss[(a, b)].append((c, d))

Start searching upper-left 4 numbers. The constraints for them are that they should be different, and should appear in the list of column-solutions.

# a b * *
# c d * *
# * * * *
# * * * *
for a, b in poss:
    for e, f in poss:
        if a == e or a == f or b == e or b == f:
            continue
        if (a, e) not in poss or (b, f) not in poss:
            continue

Next, iterate through possible values of upper-right 4 numbers. The constraints on them are that they should be all different, and should appear in the list of column-solutions.

# a b c d
# e f g h
# * * * *
# * * * *
for a, b in poss:
    for e, f in poss:
        if a == e or a == f or b == e or b == f:
            continue
        if (a, e) not in poss or (b, f) not in poss:
            continue
        for c, d in poss[(a, b)]:
            for g, h in poss[(e, f)]:
                if len(set((a, b, c, d, e, f, g, h))) != 8:
                    continue
                if (c, g) not in poss or (d, h) not in poss:
                    continue

Next, iterate through possible values of lower-left 4 numbers. The constraints on them are that they should be all different, and should appear in the list of row-solutions.

# a b c d 
# e f g h
# i j * *
# m n * *
for a, b in poss:
    for e, f in poss:
        if a == e or a == f or b == e or b == f:
            continue
        if (a, e) not in poss or (b, f) not in poss:
            continue
        for c, d in poss[(a, b)]:
            for g, h in poss[(e, f)]:
                if len(set((a, b, c, d, e, f, g, h))) != 8:
                    continue
                if (c, g) not in poss or (d, h) not in poss:
                    continue
                for i, m in poss[(a, e)]:
                    for j, n in poss[(b, f)]:
                        if len(set((a, b, e, f, i, j, m, n))) != 8:
                            continue
                        if (i, j) not in poss or (m, n) not in poss:
                            continue

Finally, iterate over the rest of the numbers. Check all constraints on them.

# a b c d 
# e f g h
# i j k l
# m n o p
for a, b in poss:
    for e, f in poss:
        if a == e or a == f or b == e or b == f:
            continue
        if (a, e) not in poss or (b, f) not in poss:
            continue
        for c, d in poss[(a, b)]:
            for g, h in poss[(e, f)]:
                if len(set((a, b, c, d, e, f, g, h))) != 8:
                    continue
                if (c, g) not in poss or (d, h) not in poss:
                    continue
                for i, m in poss[(a, e)]:
                    for j, n in poss[(b, f)]:
                        if len(set((a, b, e, f, i, j, m, n))) != 8:
                            continue
                        if (i, j) not in poss or (m, n) not in poss:
                            continue
                        for k, o in poss[(c, g)]:
                            for l, p in poss[(d, h)]:
                                solution = (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p)
                                if len(set(solution)) != 16:
                                    continue
                                if (k, l) not in poss[(i, j)] or (o, p) not in poss[(m, n)]:
                                    continue
                                # print the solution

As for larger dimensions of magic squares — I think this idea will work there too, but maybe you can find better ones.