Constructing Magic Squares of nth Powers

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, ..., k2 -- 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

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)

 xa = a11n + a12n + ... + a1rn = a21n + a22n + ... + a2rn = ... = as1n + as2n + ... + asrn

and a taxicab(n,s,r)

 xb = b11n + b12n + ... + b1sn = b21n + b22n + ... + b2sn = ... = br1n + br2n + ... + brsn

such that all possible r2s2 products aijbuv are distinct.

Let P1, ..., Ps be s Latin squares of size r × r over the elements {1,2,...,r}.

Let Q1, ..., Qr 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 Mij each of size s × r

 M11 M12 ... M1s M21 M22 ... M2s ... ... ... ... Mr1 Mr2 ... Mrs

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

(Mij)uv = (a[j, Pj[i,v]] b[i, Qi[j,u]])n (using some brackets instead of subscripts for clarity).

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

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

 Σj Σv (Mij)uv = Σj Σv (a[j, Pj[i,v]] b[i, Qi[j,u]])n (by definition) = Σj b[i, Qi[j,u]]n Σv a[j, Pj[i,v]]n (since the b[] term doesn't depend on v) = Σj b[i, Qi[j,u]]n Σv a[j,v]n (since Pj[i,v] is just a permutation of v) = xa Σj b[i, Qi[j,u]]n (by the taxicab property for a) = xa Σj b[i,j]n (since Qi[j,u] is just a permutation of j) = xaxb (by the taxicab property for b).

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

Note that the sr elements in block Mij come from products of a[j,·]n (r numbers from taxicab a) and b[i,·]n (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.

in 2 ways in 3 ways in 4 ways in 5 ways
sum of 2 squares taxicab(2,2,2) taxicab(2,2,3) taxicab(2,2,4) taxicab(2,2,5)
sum of 3 squares taxicab(2,3,2) taxicab(2,3,3) taxicab(2,3,4) taxicab(2,3,5)
sum of 4 squares taxicab(2,4,2) taxicab(2,4,3) taxicab(2,4,4) taxicab(2,4,5)
sum of 2 cubes taxicab(3,2,2) taxicab(3,2,3) taxicab(3,2,4) taxicab(3,2,5)
sum of 3 cubes taxicab(3,3,2) taxicab(3,3,3)
sum of 4 cubes taxicab(3,4,2) taxicab(3,4,3)
sum of 2 4th powers taxicab(4,2,2)

Note that a taxicab like 50 = 12 + 72 = 52 + 52 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)

 650658 = 24 + 164 + 214 + 254 = 54 + 114 + 124 + 284 = 134 + 194 + 204 + 244

and the taxicab(4,3,4)

 53809938 = 24 + 714 + 734 = 174 + 624 + 794 = 294 + 534 + 824 = 374 + 464 + 834.

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

Let P1 = P2 = P3 =

 1 2 3 4 2 3 4 1 3 4 1 2 4 1 2 3

Let Q1 = Q2 = Q3 = Q4 =

 1 2 3 2 3 1 3 1 2

Applying the method gives

 (2·2)4 (16·2)4 (21·2)4 (25·2)4 (2·71)4 (16·71)4 (21·71)4 (25·71)4 (2·73)4 (16·73)4 (21·73)4 (25·73)4
 (5·71)4 (11·71)4 (12·71)4 (28·71)4 (5·73)4 (11·73)4 (12·73)4 (28·73)4 (5·2)4 (11·2)4 (12·2)4 (28·2)4
 (13·73)4 (19·73)4 (20·73)4 (24·73)4 (13·2)4 (19·2)4 (20·2)4 (24·2)4 (13·71)4 (19·71)4 (20·71)4 (24·71)4
 (16·17)4 (21·17)4 (25·17)4 (2·17)4 (16·62)4 (21·62)4 (25·62)4 (2·62)4 (16·79)4 (21·79)4 (25·79)4 (2·79)4
 (11·62)4 (12·62)4 (28·62)4 (5·62)4 (11·79)4 (12·79)4 (28·79)4 (5·79)4 (11·17)4 (12·17)4 (28·17)4 (5·17)4
 (19·79)4 (20·79)4 (24·79)4 (13·79)4 (19·17)4 (20·17)4 (24·17)4 (13·17)4 (19·62)4 (20·62)4 (24·62)4 (13·62)4
 (21·29)4 (25·29)4 (2·29)4 (16·29)4 (21·53)4 (25·53)4 (2·53)4 (16·53)4 (21·82)4 (25·82)4 (2·82)4 (16·82)4
 (12·53)4 (28·53)4 (5·53)4 (11·53)4 (12·82)4 (28·82)4 (5·82)4 (11·82)4 (12·29)4 (28·29)4 (5·29)4 (11·29)4
 (20·82)4 (24·82)4 (13·82)4 (19·82)4 (20·29)4 (24·29)4 (13·29)4 (19·29)4 (20·53)4 (24·53)4 (13·53)4 (19·53)4
 (25·37)4 (2·37)4 (16·37)4 (21·37)4 (25·46)4 (2·46)4 (16·46)4 (21·46)4 (25·83)4 (2·83)4 (16·83)4 (21·83)4
 (28·46)4 (5·46)4 (11·46)4 (12·46)4 (28·83)4 (5·83)4 (11·83)4 (12·83)4 (28·37)4 (5·37)4 (11·37)4 (12·37)4
 (24·83)4 (13·83)4 (19·83)4 (20·83)4 (24·37)4 (13·37)4 (19·37)4 (20·37)4 (24·46)4 (13·46)4 (19·46)4 (20·46)4

which gives the 12×12 semi-magic square

 44 324 424 504 3554 7814 8524 19884 9494 13874 14604 17524 1424 11364 14914 17754 3654 8034 8764 20444 264 384 404 484 1464 11684 15334 18254 104 224 244 564 9234 13494 14204 17044 2724 3574 4254 344 6824 7444 17364 3104 15014 15804 18964 10274 9924 13024 15504 1244 8694 9484 22124 3954 3234 3404 4084 2214 12644 16594 19754 1584 1874 2044 4764 854 11784 12404 14884 8064 6094 7254 584 4644 6364 14844 2654 5834 16404 19684 10664 15584 11134 13254 1064 8484 9844 22964 4104 9024 5804 6964 3774 5514 17224 20504 1644 13124 3484 8124 1454 3194 10604 12724 6894 10074 9254 744 5924 7774 12884 2304 5064 5524 19924 10794 15774 16604 11504 924 7364 9664 23244 4154 9134 9964 8884 4814 7034 7404 20754 1664 13284 17434 10364 1854 4074 4444 11044 5984 8744 9204

The magic sum is 650658 · 53809938 = 35011866639204. The largest entry is 284·834 = 23244. 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:

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 k2, 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/(√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

pmagic = k/(2πc2m2),

or approximately 1/m2. In the 12×12 semi-magic square above, pmagic ≈ 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.

magic sum modulo 501234
correction for k=1 1.253.750.000.000.00
correction for k=2 0.311.882.810.000.00
correction for k=3 0.080.702.112.110.00
... ...............
correction for k=12 1.220.830.670.971.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

fmod = Πp prime (correction factor when expressing m as a sum of k entries modulo pe)2.

where e is the largest exponent such that pe ≤ 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, fmod ≈ 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!)2 ways of permuting the rows and columns of a k × k square. Because of symmetries, this really gives only

nperm = (k!)2 / ([k/2]! 2[k/2]+1)

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

• 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

pmagic · fmod · nperm = fmod k(k!)2 / (2πc2m2 [k/2]! 2[k/2]+1)

where

• fmod (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(m2)-time algorithm

If nperm > m2, then a solution is likely to exist. We can try random row and column permutations. After ~m2 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 nperm 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! 222) ≈ 5.64×1026 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)!2 / ((k/4)!22k/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)

 xa = a11n + a12n + ... + a1rn = a21n + a22n + ... + a2rn = ... = as1n + as2n + ... + asrn

we can

• permute the rows ai of the taxicab, and
• permute the elements aij within each row i.

It turns out that permuting the rows of the taxicabs xa and xb will permute the blocks Mij, 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 P1, ..., Ps, Q1, ..., Qr 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.

Order of the
Latin square
Number of Latin squares
with sorted first row
11
21
32
424
51344
61128960

This gives us

nlatin = A000479(r)sA000479(s)r

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

pmagic · fmod · nperm · nlatin = fmod A000479(r)sA000479(s)r rs((rs)!)2 / (2πc2m2 [rs/2]! 2[rs/2]+1).

Search algorithms

The basic algorithm is:

 For every combination of Latin squares: Search the corresponding semi-magic square for two diagonals.

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:

 Construct a list of sets of k entries with a magic sum. For every combination of Latin squares: For every set in our list: Mark the set if its entries fall on different rows and columns. If two of the marked sets obey a compatibility condition: Success!

A 16×16 magic square of 4th powers

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

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

 1950354 = 24 + 214 + 294 + 324 = 74 + 234 + 244 + 344 = 84 + 94 + 164 + 374 = 144 + 264 + 274 + 314

and

 321793923 = 34 + 854 + 974 + 1164 = 234 + 254 + 984 + 1234 = 434 + 814 + 954 + 1184 = 454 + 734 + 1064 + 1134.

One can verify with a computer that all 256 products aijbuv 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 374·1234 = 45514. 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 pmagic · fmod · nperm ≈ 6.1×10-10, which doesn't look promising. Using the Latin squares, we get a factor of 248 boost. The expected number of successes becomes pmagic · fmod · nperm · nlatin ≈ 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.

cols 1-4cols 5-8cols 9-12cols 13-16sum
rows 1-401214
rows 5-821104
rows 9-1221014
rows 13-1601124
sum444416

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

Choose(256, 16) · 0.000257 · sqrt(pmagic · fmod) ≈ 13.9 millions.

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.

 64 634 874 964 1704 17854 24654 27204 1944 20374 28134 31044 2324 24364 33644 37124
 6794 22314 23284 32984 214 694 724 1024 8124 26684 27844 39444 5954 19554 20404 28904
 9284 10444 18564 42924 7764 8734 15524 35894 6804 7654 13604 31454 244 274 484 1114
 11904 22104 22954 26354 16244 30164 31324 35964 424 784 814 934 13584 25224 26194 30074
 7364 6674 4834 464 8004 7254 5254 504 31364 28424 20584 1964 39364 35674 25834 2464
 41824 8614 28294 29524 33324 6864 22544 23524 7824 1614 5294 5524 8504 1754 5754 6004
 4004 2004 9254 2254 3684 1844 8514 2074 19684 9844 45514 11074 15684 7844 36264 8824
 30384 26464 13724 25484 38134 33214 17224 31984 7754 6754 3504 6504 7134 6214 3224 5984
 9034 864 13764 12474 17014 1624 25924 23494 19954 1904 30404 27554 24784 2364 37764 34224
 21854 22804 32304 6654 9894 10324 14624 3014 27144 28324 40124 8264 18634 19444 27544 5674
 7294 29974 6484 12964 10624 43664 9444 18884 3874 15914 3444 6884 8554 35154 7604 15204
 30684 16524 36584 31864 24704 13304 29454 25654 21064 11344 25114 21874 11184 6024 13334 11614
 13054 14404 904 9454 21174 23364 1464 15334 30744 33924 2124 22264 32774 36164 2264 23734
 17524 24824 5114 16794 10804 15304 3154 10354 27124 38424 7914 25994 25444 36044 7424 24384
 41814 18084 10174 9044 39224 16964 9544 8484 27014 11684 6574 5844 16654 7204 4054 3604
 28624 32864 27564 14844 30514 35034 29384 15824 12154 13954 11704 6304 19714 22634 18984 10224

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.

 32984 9284 10444 22954 634 964 22314 6794 22104 874 42924 26354 64 23284 18564 11904 1024 7764 8734 31324 17854 27204 694 214 30164 24654 35894 35964 1704 724 15524 16244 39444 6804 7654 814 20374 31044 26684 8124 784 28134 31454 934 1944 27844 13604 424 28904 244 274 26194 24364 37124 19554 5954 25224 33644 1114 30074 2324 20404 484 13584 10354 39224 16964 29384 23364 15334 15304 10804 35034 1464 8484 15824 21174 3154 9544 30514 23524 3684 1844 17224 7254 504 6864 33324 33214 5254 2074 31984 8004 22544 8514 38134 5524 19684 9844 3504 28424 1964 1614 7824 6754 20584 11074 6504 31364 5294 45514 7754 25994 27014 11684 11704 33924 22264 38424 27124 13954 2124 5844 6304 30744 7914 6574 12154 3014 10624 43664 29454 1624 23494 10324 9894 13304 25924 18884 25654 17014 14624 9444 24704 5674 8554 35154 13334 2364 34224 19444 18634 6024 37764 15204 11614 24784 27544 7604 11184 6004 15684 7844 3224 35674 2464 1754 8504 6214 25834 8824 5984 39364 5754 36264 7134 24384 16654 7204 18984 36164 23734 36044 25444 22634 2264 3604 10224 32774 7424 4054 19714 8264 3874 15914 25114 1904 27554 28324 27144 11344 30404 6884 21874 19954 40124 3444 21064 29524 4004 2004 13724 6674 464 8614 41824 26464 4834 2254 25484 7364 28294 9254 30384 16794 41814 18084 27564 14404 9454 24824 17524 32864 904 9044 14844 13054 5114 10174 28624 6654 7294 29974 36584 864 12474 22804 21854 16524 13764 12964 31864 9034 32304 6484 30684

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:

 24 104 244 324 584 1094 5454 13084 17444 31614 1114 5554 13324 17764 32194
 7634 10904 20714 22894 28344 7774 11104 21094 23314 28864 144 204 384 424 524
 8884 12214 19984 25534 27754 164 224 364 464 504 8724 11994 19624 25074 27254
 1054 2524 3364 6094 214 4904 11764 15684 28424 984 5954 14284 19044 34514 1194
 9804 18624 20584 25484 6864 11904 22614 24994 30944 8334 2104 3994 4414 5464 1474
 13094 21424 27374 29754 9524 2314 3784 4834 5254 1684 10784 17644 22544 24504 7844
 3244 4324 7834 274 1354 11284 15044 27264 944 4704 14524 19364 35094 1214 6054
 17864 19744 24444 6584 9404 22994 25414 31464 8474 12104 5134 5674 7024 1894 2704
 21784 27834 30254 9684 13314 4864 6214 6754 2164 2974 16924 21624 23504 7524 10344
 5444 9864 344 1704 4084 14244 25814 894 4454 10684 19684 35674 1234 6154 14764
 18694 23144 6234 8904 16914 25834 31984 8614 12304 23374 7144 8844 2384 3404 6464
 28294 30754 9844 13534 22144 7824 8504 2724 3744 6124 20474 22254 7124 9794 16024
 17694 614 3054 7324 9764 19144 664 3304 7924 10564 36834 1274 6354 15244 20324
 17164 4624 6604 12544 13864 33024 8894 12704 24134 26674 15864 4274 6104 11594 12814
 31754 10164 13974 22864 29214 15254 4884 6714 10984 14034 16504 5284 7264 11884 15184

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.

 24 104 244 324 584 1094 5454 13084 17444 31614 1114 5554 13324 17764 32194
 7634 10904 20714 22894 28344 7774 11104 21094 23314 28864 144 204 384 424 524
 8884 12214 19984 25534 27754 164 224 364 464 504 8724 11994 19624 25074 27254
 1054 2524 3364 6094 214 4904 11764 15684 28424 984 5954 14284 19044 34514 1194
 9804 6864 20584 25484 18624 11904 8334 24994 30944 22614 2104 1474 4414 5464 3994
 13094 27374 29754 21424 9524 2314 4834 5254 3784 1684 10784 22544 24504 17644 7844
 7834 4324 274 3244 1354 27264 15044 944 11284 4704 35094 19364 1214 14524 6054
 25414 22994 31464 12104 8474 5674 5134 7024 2704 1894 19744 17864 24444 9404 6584
 16924 23504 10344 7524 21624 21784 30254 13314 9684 27834 4864 6754 2974 2164 6214
 5444 344 9864 1704 4084 14244 894 25814 4454 10684 19684 1234 35674 6154 14764
 16914 23144 8904 6234 18694 23374 31984 12304 8614 25834 6464 8844 3404 2384 7144
 30754 9844 28294 13534 22144 8504 2724 7824 3744 6124 22254 7124 20474 9794 16024
 7324 17694 3054 614 9764 7924 19144 3304 664 10564 15244 36834 6354 1274 20324
 17164 13864 4624 12544 6604 33024 26674 8894 24134 12704 15864 12814 4274 11594 6104
 29214 22864 10164 31754 13974 14034 10984 4884 15254 6714 15184 11884 5284 16504 7264

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

 14244 894 25814 3744 2724 23374 31984 8504 7824 12304 6124 10684 25834 4454 8614 35094 19364 1214 2164 6754 19744 17864 4864 2974 24444 6214 6054 6584 14524 9404 1094 5454 13084 464 224 7774 11104 164 364 21094 504 31614 28864 17444 23314 5444 344 9864 13534 9844 16914 23144 30754 28294 8904 22144 4084 18694 1704 6234 24 104 244 25534 12214 7634 10904 8884 19984 20714 27754 584 28344 324 22894 7324 17694 3054 31754 22864 17164 13864 29214 10164 4624 13974 9764 6604 614 12544 7924 19144 3304 15254 10984 33024 26674 14034 4884 8894 6714 10564 12704 664 24134 19684 1234 35674 9794 7124 6464 8844 22254 20474 3404 16024 14764 7144 6154 2384 4904 11764 15684 3784 4834 11904 8334 2314 5254 24994 1684 984 22614 28424 30944 7834 4324 274 7524 23504 25414 22994 16924 10344 31464 21624 1354 8474 3244 12104 5954 14284 19044 17644 22544 2104 1474 10784 24504 4414 7844 1194 3994 34514 5464 27264 15044 944 9684 30254 5674 5134 21784 13314 7024 27834 4704 1894 11284 2704 1054 2524 3364 21424 27374 9804 6864 13094 29754 20584 9524 214 18624 6094 25484 1114 5554 13324 25074 11994 144 204 8724 19624 384 27254 32194 524 17764 424 15244 36834 6354 16504 11884 15864 12814 15184 5284 4274 7264 20324 6104 1274 11594

The magic sum is 794179 · 292965218 = 232666823866022. The largest entry is 294·1274 = 36834. 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.

powersemi-magic squaremagic sumcfmodpmagic · fmod · npermpmagic · fmod · nperm · nlatinmagic square
412×12 FL3.50×10131.8725.22.8×10-146.2×10-9none (FL)
15×15 CB2.33×10142.0017.72.6×10-1020found (FL)
16×16 FL6.28×10141.7814.16.1×10-1067found (FL)
525×25 FL4.75×10252.321.222.5×10-144.7×1017?
25×25 CB8.28×10252.653.021.5×10-142.9×1017?

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!

page last updated: November 19, 2011