Constructing Magic Squares of nth Powers

This page describes the methods I used to construct:

• a 12×12 semi-magic square of 4th powers, the first known of exactly that size - April 2010,
• a 16×16 magic square of 4th powers, the smallest known at the time - December 2010,
• a 15×15 magic square of 4th powers, the new smallest known, beating my own record - September 2011.

Contents

• Definitions
• Background
• The taxicab method (semi-magic version)
  • Some taxicabs
  • A 12×12 semi-magic square of 4th powers
  • New semi-magic squares from other people
• Obtaining magic diagonals given a semi-magic square
  • A random model
  • Modular correction
  • Permuting rows and columns
  • An O(m²)-time algorithm
  • An O(m)-time algorithm
  • An O(√m)-time, O(√m)-space algorithm
• Obtaining magic diagonals given compatible taxicabs
  • Rearranging the taxicabs
  • Choosing Latin squares
  • Search algorithms
  • A 16×16 magic square of 4th powers
  • A 15×15 magic square of 4th powers
  • Other candidates for magic diagonals
• Further generalization

Definitions

• A magic square is a square grid of distinct positive integers such that each row, column, and the two diagonals have the same sum. This is the definition used for the Enigmas on Magic Squares challenge. Note that definitions can vary: sometimes negative integers are allowed; sometimes the integers are required to be the numbers 1, 2, 3, ..., k² -- but not here. See Zaslav's comment on the Wikipedia magic square talk page about this.
• A semi-magic square is a magic square but with no requirement on the diagonal sums.
• A generalized taxicab(n,i,j) number is an integer x which can be expressed as the sum of i positive nth-powers in j different ways. Wikipedia defines it as the minimum such x, but here I allow any solution.

Background

In 2006, Lee Morgenstern found methods to construct

• a 4×4 semi-magic square of nth powers by using two taxicab(n,2,2) numbers.
• a 6×6 semi-magic square of nth powers by using a taxicab(n,3,2) number and a taxicab(n,2,3) number.
• a 9×9 semi-magic square of nth powers by using two taxicab(n,3,3) numbers.

To me these are not individual results, but rather special cases of a more general method to construct a (rs)×(rs) semi-magic square by using a taxicab(n,r,s) and a taxicab(n,s,r). My first goal is to explain the generalization, which I call the taxicab method.

The taxicab method (semi-magic version)

Choose a taxicab(n,r,s)

and a taxicab(n,s,r)

such that all possible r²s² products a_ijb_uv are distinct.

Let P₁, ..., P_s be s Latin squares of size r × r over the elements {1,2,...,r}.

Let Q₁, ..., Q_r be r Latin squares of size s × s over the elements {1,2,...,s}.

Construct a (rs)×(rs) square M consisting of r × s blocks M_ij each of size s × r

where element (u,v) of block M_ij is given by

(M_ij)_uv = (a[j, P_j[i,v]] b[i, Q_i[j,u]])ⁿ (using some brackets instead of subscripts for clarity).

Then M is a (rs)×(rs) semi-magic square of nth powers, with magic sum x_ax_b.

Proof.
To compute the sum of every element in a row, we fix i and u, and sum over j and v:

So every row of the square M sums to x_ax_b. The proof that every column sums to x_ax_b is similar. End of proof.

Note that the sr elements in block M_ij come from products of a[j,·]ⁿ (r numbers from taxicab a) and b[i,·]ⁿ (s numbers from taxicab b). The specific arrangement of those sr numbers is determined by the Latin squares.

Some taxicabs

Below are links to The On-Line Encyclopedia of Integer Sequences corresponding to some generalized taxicabs. The sequences only give the taxicab sums, not how to realize the sums, but you can find that information in at least one reference when the link is highlighted.

Note that a taxicab like 50 = 1² + 7² = 5² + 5² can't be used to create a semi-magic square because its components aren't distinct. We can also discard a taxicab when the gcd of its components isn't 1 because the gcd will divide the semi-magic square entries.

A 12×12 semi-magic square of 4th powers

Let n = 4, r = 4, and s = 3.

I choose the taxicab(4,4,3)

and the taxicab(4,3,4)

One can verify with a computer that all 144 products a_ijb_uv are distinct.

Let P₁ = P₂ = P₃ =

Let Q₁ = Q₂ = Q₃ = Q₄ =

Applying the method gives

which gives the 12×12 semi-magic square

The magic sum is 650658 · 53809938 = 35011866639204. The largest entry is 28⁴·83⁴ = 2324⁴. This is the smallest possible magic sum and smallest possible largest entry using compatible taxicab(4,4,3) and taxicab(4,3,4). The diagonal sums are 12005752179109 and 45297006668644 which are not magic.

New semi-magic squares from other people

After I published my taxicab method formula on this webpage, Christian Boyer created:

• a 25×25 semi-magic square of 5th powers - March 2011,
• a 15×15 semi-magic square of cubes with the smallest possible magic sum using taxicabs - May 2011,
• a 15×15 semi-magic square of 4th powers with the smallest possible magic sum using taxicabs - May 2011.

Obtaining magic diagonals given a semi-magic square

A random model

Semi-magic squares produced by the taxicab method are unlikely to be fully magic. For a k × k semi-magic square with magic sum m, the entries can be thought as random with mean m/k and standard deviation c·m/k where in practice c is close to 1. For example, c ≈ 0.577 if the entries are the numbers 1 through k², and c ≈ 1.867 for the 12×12 semi-magic square above. Under this random model, diagonal sums have mean m and standard deviation c·m/√k. Assuming a normal distribution, the probability that a given diagonal sum equals m is √k/(√2πc·m). Typically this expression is dominated by m, so for simplicity I will say that the probability is ~1/m. To get two correct diagonal sums, we need to square this probability, which gives

or approximately 1/m². In the 12×12 semi-magic square above, p_magic ≈ 4.5×10^-28.

Modular correction

The random model above requires a correction because magic square entries are not always uniformly distributed modulo prime powers, so a sum of k entries is not necessarily uniformly distributed modulo prime powers either. The effect can be surprisingly important. I define the correction factor modulo M to be M times the probability of obtaining the magic sum m as a sum of k random entries modulo M.

As an example, the table below gives correction factors modulo 5 assuming a square with entries chosen from the 12×12 semi-magic square above. The actual square has k=12 and m≡4 (mod 5), so its correction factor modulo 5 is 1.30.

The complete correction factor takes into accounts all prime powers, and must be squared because we're matching 2 diagonals. It can be defined as

where e is the largest exponent such that p^e ≤ 1024 (a bound to make sure that the effect is local). The modular correction tells how much more likely it is to get 2 diagonals summing to m because of modular considerations. For the 12×12 semi-magic square above, f_mod ≈ 25.2, primarily due to effects modulo 5 and modulo 16.

Permuting rows and columns

Given a semi-magic square, we can generate more semi-magic squares by permuting its rows and columns any way we like. This gives us more than one chance of getting correct diagonals. There are (k!)² ways of permuting the rows and columns of a k × k square. Because of symmetries, this really gives only

independent chances of success (for k ≥ 2), where [] denotes the floor function. For k=12, this is 4.98×10¹².

• The factor [k/2]! comes from the fact that the first [k/2] rows can be permuted without changing the diagonal sums, as long as the last [k/2] rows are permuted symmetrically, and the columns are permuted in the same way as the rows.
• The factor 2^[k/2] comes from the fact that we can swap rows x and k+1-x without changing the diagonal sums, as long as we swap the same columns.
• The extra factor 2 comes from the fact that we can interchange the two diagonals with a reflection (for k ≥ 2).

Combining the last three formulas, the expected number of magic diagonal pairs that can be obtained (counting equivalent cases only once) is

where

• f_mod (typically in the range 1 - 25) is the modular correction factor of the square
• k is the order of the square
• c (typically in the range 1.7 - 2.7) is the coefficient of variation of the entries
• m is the magic sum.

An O(m²)-time algorithm

If n_perm > m², then a solution is likely to exist. We can try random row and column permutations. After ~m² attempts we should find a solution.

An O(m)-time algorithm

Note that if we apply the same permutation to both rows and columns, then elements from the main diagonal stay in the main diagonal, and so their sum stays unchanged. This suggests the following algorithm:

1. Apply random permutations of rows (or columns) until we get a correct main diagonal. This should take ~m attempts.
2. Generate a random permutation and apply it to both rows and columns until we get a correct anti-diagonal. This should take ~m attempts.

Mathematically, this is splitting n_perm as k!/2 × k!/([k/2]! 2^[k/2]). We need k!/([k/2]! 2^[k/2]) > m for a solution to step 2 to be likely to exist.

This is apparently the method that was used by Hugo Pfoertner to find magic diagonals in a 44×44 semi-magic square of 4th powers in October 2007. In his case, m=123784061778806 and k=44, so 44!/(22! 2²²) ≈ 5.64×10²⁶ which is plenty large for the algorithm to work.

An O(√m)-time, O(√m)-space algorithm

If we have plenty of permutations available, then we can add constraints to the permutations and still be ok. The trick is to split the indices {1,2,...,k} into two sets S and T. Generate ~√m permutations S→S and store them keyed by the sum of diagonal entries restricted to S. Similarly, generate ~√m permutations T→T and store them keyed by m - (sum of diagonal entries restricted to T). Look for a collision between the two sets to obtain a diagonal with sum m.

Doing the anti-diagonal is tricky as each set S and T should be symmetric (x is in S iff k+1-x is in S). One can choose S = {1,2,...,q} union {k+1-q, ..., k-1, k} where q ≈ k/4, and let T be the remaining middle part. The condition for this algorithm to work is roughly (k/2)!² / ((k/4)!²2^k/2) > m.

Obtaining magic diagonals given compatible taxicabs

In the previous section I discussed turning a given semi-magic square into a fully magic square by permuting rows and columns. But if the semi-magic square comes from the taxicab method, there is even more flexibility.

Rearranging the taxicabs

Note that in the taxicab(n,r,s)

we can

• permute the rows a_i of the taxicab, and
• permute the elements a_ij within each row i.

It turns out that permuting the rows of the taxicabs x_a and x_b will permute the blocks M_ij, which doesn't give us anything new since we're already planning to permute the rows and columns of M individually. Permuting the elements within each row of the taxicabs does give us something new, but the same permutation can be accomplished by choosing Latin squares, so in the end there is no need to rearrange the taxicabs.

Choosing Latin squares

We can assume that the first row of each Latin square P₁, ..., P_s, Q₁, ..., Q_r is ordered, because any other ordering can be obtained later by permuting the rows and columns of M. We do not assume that the first column of each Latin square is ordered because a similar claim doesn't hold for columns. We're interested in the number of such Latin squares, which is given by the integer sequence A000479.

This gives us

possible combinations. For a given compatible taxicab pair taxicab(n,r,s) and taxicab(n,s,r), the expected number of possible magic diagonal pairs becomes

Search algorithms

The basic algorithm is:

This works fine for small squares, but for large squares it is too expensive to do O(√m) work for each combination of Latin squares. The trick is to do some precomputation in order to speed up the inner loop. The revised algorithm is:

A 16×16 magic square of 4th powers

Let n = 4, r = 4, and s = 4.

I choose the two taxicabs(4,4,4)

and

One can verify with a computer that all 256 products a_ijb_uv are distinct. As we know by now, this gives a 16×16 semi-magic square with magic sum m = 1950354 · 321793923 = 627612064898742. The largest entry is 37⁴·123⁴ = 4551⁴. This is the smallest possible magic sum and smallest possible largest entry using two compatible taxicabs(4,4,4).

Using the formulas above, the expected number of successes if we search all possible permutations is p_magic · f_mod · n_perm ≈ 6.1×10^-10, which doesn't look promising. Using the Latin squares, we get a factor of 24⁸ boost. The expected number of successes becomes p_magic · f_mod · n_perm · n_latin ≈ 67, which is enough.

The numbers were tight and I could not afford to add constraints on the permutations.

The first step of the algorithm was to search for sets of 16 entries such that

1. their sum equals the magic sum, and
2. their positions can potentially form a diagonal after rearrangement.

To explain what I mean by (2), recall from the taxicab method (semi-magic version) that the 256 entries are arranged in 16 blocks of 16 numbers, and that the Latin squares shuffle numbers within blocks and never swap numbers across blocks. Consider a diagonal in a 16×16 square. Note that if we permute the rows and columns of the square, then the diagonal entries get shuffled, but we still get 1 diagonal entry per row, and 1 diagonal entry per column. Now if we further shuffle the diagonal entries using Latin squares, then those properties don't necessarily hold, but what remains true is that the top 4 rows together contain 4 diagonal entries, and similarly for the next groups of 4 rows, and same for columns.

To illustrate this, the table below gives the number of diagonal entries that fall into each block of 4×4 entries for some arbitrary permitted shuffling. Those counts can be anything from 0 to 4, but note that the sums across rows and columns are always 4.

It turns out that among all Choose(256, 16) possible ways to choose a diagonal's entries, only about 0.0257% of the choices satisfy those sum conditions. I searched for those only, which reduced the work considerably.

The expected number of sets of 16 entries that satisfy both (1) and (2) is

I used a collision method to search for about 67% of the possibilities. I found 9429609 sets, which is pretty close to what the random model predicted.

Next, for each combination of Latin squares that I tried, there were typically ~850 sets that fell on distinct rows, and an expected 0.08 sets that also fell on distinct columns. Each time there was a small probability of finding two compatible magic diagonals in the same square. Eventually I found an example.

By looking at the pattern of colors, can you guess what the compatibility condition is? This allowed me to permute the rows and columns to get a fully magic square.

At the time of discovery, this 16×16 square was the smallest known magic square of 4th powers. The previous record was the 36×36 pentamagic square of Li Wen, found in 2008.

A 15×15 magic square of 4th powers

Starting from Christian Boyer's 15×15 semi-magic square of 4th powers, I applied the same method above to find magic diagonals.

Here's Christian's published semi-magic square, with the elements that will become diagonals highlighted:

Here's the square after the application of the Latin squares that I selected. This demontrates the importance of Latin squares when searching for diagonals.

Here's the fully magic square, after permuting the rows and columns:

The magic sum is 794179 · 292965218 = 232666823866022. The largest entry is 29⁴·127⁴ = 3683⁴. This is the smallest possible magic sum and smallest possible largest entry using compatible taxicab(4,5,3) and taxicab(4,3,5). At the time of writing, this 15×15 square is the smallest known magic square of 4th powers. The previous record was my 16×16 square described earlier.

Other candidates for magic diagonals

Here are semi-magic squares that could lead to new records if their entries can be rearranged to form 2 magic diagonals:

• Christian Boyer's 25x25 semi-magic square of 5th powers,
• the 25×25 semi-magic square of 5th powers with the smallest possible magic sum, using the two taxicabs(5,5,5)

  1189260043200	=	225 + 585 + 745 + 1615 + 2555
  	=	285 + 625 + 1555 + 1655 + 2505
  	=	435 + 1495 + 1725 + 1765 + 2405
  	=	795 + 1195 + 1865 + 1955 + 2315
  	=	1235 + 1345 + 1875 + 2055 + 2215

  and

  39931433836700	=	105 + 2235 + 2925 + 3825 + 4935
  	=	335 + 1915 + 2995 + 3855 + 4925
  	=	355 + 1995 + 2715 + 4125 + 4835
  	=	1135 + 2275 + 3655 + 4165 + 4595
  	=	2635 + 3325 + 3895 + 4065 + 4305.

Below is a calculation of the expected number of essentially distinct squares with magic diagonals according to the random model for these two candidates and for the three semi-magic squares I discussed earlier, sorted by power and by size.

The results show that a 25×25 magic square of 5th powers can almost certainly be obtained from the respective semi-magic squares, but to succeed we almost certainly have to use the flexibility afforded by the Latin squares.

Further generalization

In 2006, Christian Boyer found a method to construct an 8×8 semi-magic square of nth powers by using three taxicab(n,2,2) numbers. This hints at a general method to construct a (rst)×(rst) semi-magic square of nth powers by using a taxicab(n,r,s), a taxicab(n,s,t), and a taxicab(n,t,r). To push the idea further, maybe if we write k as the product of d factors, it is possible to create a semi-magic square of size k × k using d taxicabs? What about the Latin squares, how do they interact with this? Someone can try to figure it out!