Numerous but Rare: An Exploration of Magic Squares

Introduction

Making magic squares is a popular form of mathematical recreation. It is also used in classrooms as an elementary mathematical exercise. The classic (or ordinary) magic square of order n is defined as follows: Placing the numbers 1, 2, ⋯n ² in a square array using each number once, if all the sums of the numbers in each row, column and diagonal give the same value, $M_n=\frac{1}{2}n(n^2+1)$, the array makes a magic square. M _n is called the magic number. Besides the classic magic squares, there are many variations, and some rigorous results have been found for them. But not much is known about classic magic squares [1]. In this paper, we focus on classic magic squares.

There are some algorithms for making magic squares of any size. They, however, provide some special classes of magic squares, which gives rise to the question: Among all the possible arrangements of numbers in a square of a given size, how many of them form magic squares? Putting the question in another form: Is there any chance of making a magic square by putting numbers randomly in a square? It may be surprising to know that the exact number of possible magic squares is so far only known up to order 5. There is currently no hope of exact enumeration for a larger system. In this paper, we apply a Monte Carlo method to this problem, and estimate the number of the magic squares of each size up to order 30.

To state the problem explicitly, we consider classic magic squares order n. Let N _n denote the total number of magic squares. Since possible configurations increase as n ²!, counting the magic squares rapidly becomes more difficult. Currently, only the following three values of N _n are known exactly: N ₃ = 1, N ₄ = 880, and N ₅ = 275,305,224 [1], where the eight equivalent magic squares that can be transformed into each other by rotation and reflection are counted as one.

For larger squares, we need to employ statistical approaches for estimating the number of magic squares. There have been two studies in this direction. Pinn and Wieczerkowski applied the exchange Monte Carlo method (EMC) [2] to this problem [3] and estimated N ₆ and N ₇. Their results are N ₆ = 1.7745(16) × 10¹⁹ and N ₇ = 3.760(52) × 10³⁴, where the digits in parentheses indicate the statistical error of the lowest digits. Trump proposed a more efficient method, called Monte Carlo Backtracking, and estimated N _n for n ≤ 20. [4].

EMC belongs to a family of extended ensemble Monte Carlo methods [5]. Extended ensemble Monte Carlo methods were initially developed in the field of statistical physics, and have found a wide field of applications beyond their original scope. They are especially suitable for estimating the probability of occurrence of very rare events. The work by Pinn and Wieczerkowski is one of the earliest applications of the extended ensemble Monte Carlo methods outside the field of physics. In this paper, we use the Multicanonical Monte Carlo method (MMC) [6], which also belongs to a family of extended ensemble methods. There have also been some studies that used EMC for counting solutions for mathematical puzzles such as the N-Queen problem [7][8]. But the MMC has not been used for problems of this type. Compared to the EMC, which requires a trick for counting the number of configurations that satisfy some specific conditions, the MMC provides the estimates of the number straightforwardly. We thus consider the MMC to be more suitable for problems of this sort than the EMC.

Methods

Let us consider magic squares of order n. In order to apply the MMC, we define the energy E(C) of a configuration of numbers C as follows:

where S _i, S _j, and S _k are the sums of the numbers for the ith column, that for the jth row, and that for the kth diagonal. E(C) is zero if and only if C is a magic square, and it takes a positive value otherwise. Thus, the lowest energy (E = 0) configurations provide magic squares. In other words, finding problem of magic squares are formulated as a combinatorial optimization problem. The number of optimal configurations are very large in this case, and we estimate the number using the MMC.

The MMC was proposed by Berg and Neuhaus in the field of statistical physics to overcome slow convergence of Metropolis-type Markov chain Monte Carlo methods when applied to the sampling of low temperature states of complex systems [6]. In contrast to the Metropolis method which provides the canonical ensemble, namely the ensemble of configurations at the thermal equilibrium of a given temperature T as the steady-state distribution of the Markov process, the MMC aims to give a flat energy distribution over the entire energy range. This flatness enables us to estimate the number of configurations of any energy. For that purpose, a predetermined weight function W(E) is used in the MMC instead of the canonical weight e ^−E/T used in the Metropolis method. W(E) is prepared so that the energy histogram H(E) obtained by Monte Carlo sampling is sufficiently flat. Since H(E) is proportional to the product of W(E) and the number of states g(E) having energy E, we can then estimate relative values of g(E) from W(E) and H(E) as

The appearance probability of energy E = ɛ in randomly arranged configurations of numbers from 1 to n ² is estimated by

where the summation in the denominator is taken over all the possible energies. Thus the appearance probability of magic square P _n is given by P _n = P(E = 0). Since the total number of configurations is n ²!, we can also estimate the total number of magic squares N _n by P _n × n ²!/8. It should be noted that, in principle, MMC gives statistically unbiased estimates for the number of magic squares.

To determine W(E), we carried out preliminary runs using the Wang–Landau method [9], in which W(E) is updated at each Monte Carlo trial until it finally gives a sufficiently flat histogram H(E). We then fixed W(E) and carried out long measurement runs using the entropic sampling method [10], which is equivalent to MMC in the present study because we assigned a single value of energy to each bin of the histogram. We made independent measurement runs many times for each n, and took averages of N _n and $N_n^2$ over them. The statistical error of N _n was then estimated as three times the standard error:

where ⟨⋅⟩ is the mean value of ⋅ and t is the number of the measurement runs.

Only the sequential transposition of adjacent numbers, 1 with 2, then 2 with 3, ⋯, n ²−1 with n ² were used as an elementary process of the Monte Carlo trial, following Ref. [3]. By this method, we can avoid a large energy difference between successive configurations in Markov chains, which causes inefficiencies in Monte Carlo methods. We employed Mersenne-Twister as the pseudo-random-number generator [11].

The number of measurement runs and the length of each run are different for each n. For the largest system with n = 30, for example, we made 40 independent measurement runs of 1.1 × 10¹² Monte Carlo trials each. Flatness of the histogram was confirmed by a long independent run that was four times longer than the measurement run. The number of configurations in each bin of the histogram falls within the range from 0.93 to 1.01 of the mean value, which we decided as sufficiently flat.

Results

Fig 1 shows a semi-log plot of the n dependence of P _n. Our estimates of P _n and N _n up to n = 30 are listed in Table 1. Exact values for of N _n for n = 3,4 and 5 and the previous estimates of N _n by Trump for 6 ≤ n ≤ 20 are also shown for comparison. We obtained N ₃₀ = 6.56(29) × 10²⁰⁵⁶. In contrast to this extremely large value, its appearance probability, P ₃₀, is smaller than 10⁻²¹². Thus the magic squares are numerous but rare.

The largest size of n = 30 is much larger than that which has been calculated previously. The estimates of N ₃, N ₄, and N ₅ agree with the exact values within the statistical error, and the estimates up to n = 17 are consistent with Trump’s values within the statistical error. However, there are appreciable discrepancies between the present results and those of Trump for n = 18 and 19; our values are larger for these sizes. In fact, Trump himself pointed out that the true values for N ₁₈ and N ₁₉ might be smaller than his estimates based on his own extrapolation formula. We thus think that our estimates are reliable. He also gave estimates of N _n for n > 20 obtained by the abovementioned extrapolation formula. Compared to the present estimates, his extrapolation results have two-digits accuracy up to n = 30.

As seen in Fig 1, the appearance probability of magic squares P _n decreases rapidly with n. In other words, magic squares become rarer rapidly as n increases. This raises the question: how fast does their appearance probability decrease? It clearly decreases faster than the exponential function exp(−an) with constant a. On the other hand, since the number of possible configurations is n ²!, P _n should decrease slower than exp(−n ²logn). Considering these facts, we first tried to fit logP _n for n ≥ 10 by the second-order polynomial. But the reduced χ ² was larger than 6000 and thus the fitting function was inappropriate. Next we tried functions including logn. The fitting function (An + B)logn + Cn + D with the constants A, B, C and D gave A = −4.99 ± 0.07 with the reduced χ ² = 2.55 and the fitting function (En + F)log(n + G) + H with the constants E, F, G and H gave E = −4.880 ± 0.008 with the reduced χ ² = 2.42. Fitted curves using these two functions are shown in Fig 1. The two curves are virtually indistinguishable on this scale except for very small values of n. Other functions we tried gave larger values of reduced χ ². The reduced χ ² for both functions, however, are still large, and the fittings are not fully satisfactory. We consider the reason is that n = 10 is not yet at the asymptotic region. In fact, when we tried to fit the same functions to data only for larger sizes, n ≥ 20, we obtained the reduced χ ² = 1.30 and 1.37, respectively. Although the errors of the parameters are large, A = −4.6 ± 0.5 and E = −4.88 ± 0.06, reduced χ ² for both functions are rather satisfactory.

Fig 2 shows the ratio P _n/exp{(An + B)logn + Cn + D} and P _n/exp{(En + F)log(n + G) + H}. Both functions seem to express P _n equally well. In any case, since the slope of logP _n varies slowly, it is difficult to determine the appropriate functional form from the present results. We need P _n for much larger systems to provide a solid conclusion. Even so, we may conjecture that the logP _n decreases asymptotically as an log n with a ≃ 5.

In this paper, we applied the Multicanonical Monte Carlo method to a combinatorial optimization problem by defining appropriate an energy function E(C). The MMC directly gives the number of the optimal configurations from the histogram of the lowest energy configurations. The present work demonstrates that the MMC is a powerful tool for counting rare configurations of combinatorial problems. We can estimate the appearance probabilities of the optimal configurations as small as 10⁻²¹².

Data Availability

All relevant data are within the paper.

Funding Statement

AK is supported by Global COE Program (Core Research and Engineering of Advanced Materials Interdisciplinary Education Center for Materials Science), MEXT, Japan (http://www.gcoe.mp.es.osaka-u.ac.jp).