Magic squares are a fascinating mathematical challenge that has intrigued mathematicians for centuries. Given a positive (and possibly large) integer \( n \) , one of the main challenges that still remains is to find, within a reasonable computational time, a magic square of order \( n \) , that is, a square matrix of order \( n \) with unique integers from \( a_{\min } \) to \( a_{\max } \) , such that the sum of each row, column, and diagonal equals a constant \( \mathcal {C}(A) \) . In this work, we first present an integer constraint satisfaction problem for constructing a magic square of order \( n \) . Nonetheless, the solution time of this problem grows exponentially as the order increases. To overcome this limitation, we also propose a fast approach that constructs magic squares depending on whether \( n \) is odd, singly even, or doubly even. Moreover, we provide a proof of the correctness of this novel approach. Our numerical results show that our method can construct magic squares of order up to 70 000 in less than 140 seconds, demonstrating its efficiency and scalability.
