Abstract

Citation:Kitajima A, Kikuchi M (2015) Numerous but Rare: An Exploration of Magic Squares. PLoS ONE 10(5): e0125062. https://doi.org/10.1371/journal.pone.0125062

Academic Editor:Eduardo G. Altmann, Max Planck Institute for the Physics of Complex Systems, GERMANY

Received:December 21, 2014;Accepted:March 20, 2015;Published: May 14, 2015

Copyright: © 2015 Kitajima, Kikuchi. This is an open access article distributed under the terms of theCreative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited

Data Availability:All relevant data are within the paper.

Funding: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).

Competing interests: The authors have declared that no competing interests exist.

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 ordern 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 ordern. LetN_n denote the total number of magic squares. Since possible configurations increase asn²!, counting the magic squares rapidly becomes more difficult. Currently, only the following three values ofN_n are known exactly:N₃ = 1,N₄ = 880, andN₅ = 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 estimatedN₆ andN₇. Their results areN₆ = 1.7745(16) × 10¹⁹ andN₇ = 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 estimatedN_n forn ≤ 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 ordern. In order to apply the MMC, we define theenergyE(C) of a configuration of numbersC as follows:(1)whereS_i,S_j, andS_k are the sums of the numbers for theith column, that for thejth row, and that for thekth diagonal.E(C) is zero if and only ifC 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 temperatureT 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 functionW(E) is used in the MMC instead of the canonical weighte^−E/T used in the Metropolis method.W(E) is prepared so that the energy histogramH(E) obtained by Monte Carlo sampling is sufficiently flat. SinceH(E) is proportional to the product ofW(E) and the number of statesg(E) having energyE, we can then estimate relative values ofg(E) fromW(E) andH(E) as(2)

The appearance probability of energyE =ɛ in randomly arranged configurations of numbers from 1 ton² is estimated by(3)where the summation in the denominator is taken over all the possible energies. Thus the appearance probability of magic squareP_n is given byP_n =P(E = 0). Since the total number of configurations isn²!, we can also estimate the total number of magic squaresN_n byP_n ×n²!/8. It should be noted that, in principle, MMC gives statistically unbiased estimates for the number of magic squares.

To determineW(E), we carried out preliminary runs using the Wang–Landau method [9], in whichW(E) is updated at each Monte Carlo trial until it finally gives a sufficiently flat histogramH(E). We then fixedW(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 eachn, and took averages ofN_n and$N_n^2$ over them. The statistical error ofN_n was then estimated as three times the standard error:(4)where ⟨⋅⟩ is the mean value of ⋅ andt is the number of the measurement runs.

Only the sequential transposition of adjacent numbers, 1 with 2, then 2 with 3, ⋯,n²−1 withn² 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 eachn. For the largest system withn = 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 then dependence ofP_n. Our estimates ofP_n andN_n up ton = 30 are listed inTable 1. Exact values for ofN_n forn = 3,4 and 5 and the previous estimates ofN_n by Trump for 6 ≤n ≤ 20 are also shown for comparison. We obtainedN₃₀ = 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Semi-log plot of the appearance probabilityP_n of magic squares (•).

P_n decreases faster than exponentially with the sizen. Two fitted functions are also shown: exp((An +B)ln(n) +Cn +D) (solid line) and exp((En +F)ln(n +G) +H) (dotted line) withA = −4.99 andE = −4.88. We usedP_n ofn ≥ 10 for the fitting. Enlarged plot forn < 6 is shown in the inset, in which difference of two functions are visible.

https://doi.org/10.1371/journal.pone.0125062.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Estimated number and qppearance probability of magic squares.

https://doi.org/10.1371/journal.pone.0125062.t001

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

As seen inFig 1, the appearance probability of magic squaresP_n decreases rapidly withn. In other words, magic squares become rarer rapidly asn increases. This raises the question: how fast does their appearance probability decrease? It clearly decreases faster than the exponential function exp(−an) with constanta. On the other hand, since the number of possible configurations isn²!,P_n should decrease slower than exp(−n²logn). Considering these facts, we first tried to fit logP_n forn ≥ 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 constantsA,B,C andD gaveA = −4.99 ± 0.07 with the reducedχ² = 2.55 and the fitting function (En +F)log(n +G) +H with the constantsE,F,G andH gaveE = −4.880 ± 0.008 with the reducedχ² = 2.42. Fitted curves using these two functions are shown inFig 1. The two curves are virtually indistinguishable on this scale except for very small values ofn. 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 thatn = 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 andE = −4.88 ± 0.06, reducedχ² for both functions are rather satisfactory.

Fig 2 shows the ratioP_n/exp{(An +B)logn +Cn +D} andP_n/exp{(En +F)log(n +G) +H}. Both functions seem to expressP_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 needP_n for much larger systems to provide a solid conclusion. Even so, we may conjecture that the logP_n decreases asymptotically asan logn witha ≃ 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The ratio of the appearance probabilityP_n to two fitted functions.

P_n/exp{(An +B)ln(n) +Cn +D}} (•) andP_n/exp{(En +F)ln(n +G) +H} (×) withA = −4.99 andE = −4.88. Both functions seem to expressP_n equally well.

https://doi.org/10.1371/journal.pone.0125062.g002

In this paper, we applied the Multicanonical Monte Carlo method to a combinatorial optimization problem by defining appropriate an energy functionE(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⁻²¹².

Acknowledgments

We are grateful to Koji Hukushima, Yukito Iba, and Nen Saito for valuable discussions and continuous encouragement. Numerical experiments were carried out on PC clusters at the Cybermedia Center, Osaka University.

Author Contributions

Conceived and designed the experiments: AK MK. Performed the experiments: AK. Analyzed the data: AK MK. Contributed reagents/materials/analysis tools: AK MK. Wrote the paper: AK MK.

References

1. 1.Pickover C (2001) The Zen of Magic Squares, Circles, and Stars. Princeton University Press.
2. 2.Hukushima K, Nemoto K (1996) Exchange Monte Carlo Method and Application to Spin Glass Simulations. J Phys Soc Jpn 65 (1996): 1604–1608.
   • View Article
   • Google Scholar
3. 3.Pinn K, Wieczerkowski C (1998) Number of Magic Squares from Parallel Tempering Monte Carlo. Int J Mod Phys C9: 541–546.
   • View Article
   • Google Scholar
4. 4.Trump W (2005). Available:http://www.trump.de/magic-squares. Accessed 17 December 2014.
5. 5.Iba Y (2001) Extended ensemble Monte Carlo. Int J Mod Phys C12: 623–656.
   • View Article
   • Google Scholar
6. 6.Berg BA, Neuhaus T (1991) Multicanonical algorithms for first order phase transitions. Phys Lett B267: 249–253.
   • View Article
   • Google Scholar
7. 7.Hukushima K (2002) Extended ensemble Monte Carlo approach to hardly relaxing problems. Computer physics communications 147: 77–82.
   • View Article
   • Google Scholar
8. 8.Zhang C, Ma J (2009) Counting solutions for the N-queens and Latin-square problems by Monte Carlo simulations. Phys Rev E79: 016703.
   • View Article
   • Google Scholar
9. 9.Wang F, Landau DP (2001) Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States. Phys Rev Lett 86: 2050–2053. pmid:11289852
   • View Article
   • PubMed/NCBI
   • Google Scholar
10. 10.Lee J (1993) New Monte Carlo Algorithm: Entropic Sampling. Phys Rev Lett 71: 211–214. pmid:10054892
    • View Article
    • PubMed/NCBI
    • Google Scholar
11. 11.Matsumoto M, Nishimura T (1998) Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans on Modeling and Computer Simulation 8: 3–30.
    • View Article
    • Google Scholar