Constructing Magic Squares of nth Powers
This page describes the methods I used to construct:
Contents
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^{2}  but not here. See Zaslav's comment on the Wikipedia magic square talk page about this.
 A semimagic 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 nthpowers 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) semimagic 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 (semimagic version)
Choose a taxicab(n,r,s)
x_{a}  =  a_{11}^{n} + a_{12}^{n} + ... + a_{1r}^{n} 
 =  a_{21}^{n} + a_{22}^{n} + ... + a_{2r}^{n} 
 =  ... 
 =  a_{s1}^{n} + a_{s2}^{n} + ... + a_{sr}^{n} 
and a taxicab(n,s,r)
x_{b}  =  b_{11}^{n} + b_{12}^{n} + ... + b_{1s}^{n} 
 =  b_{21}^{n} + b_{22}^{n} + ... + b_{2s}^{n} 
 =  ... 
 =  b_{r1}^{n} + b_{r2}^{n} + ... + b_{rs}^{n} 
such that all possible r^{2}s^{2} products a_{ij}b_{uv} are distinct.
Let P_{1}, ..., P_{s} be s Latin squares of size r × r over the elements {1,2,...,r}.
Let Q_{1}, ..., 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
M_{11}  M_{12}  ...  M_{1s} 
M_{21}  M_{22}  ...  M_{2s} 
...  ...  ...  ... 
M_{r1}  M_{r2}  ...  M_{rs} 
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]])^{n}
(using some brackets instead of subscripts for clarity).
Then M is a (rs)×(rs) semimagic square of nth powers, with magic sum x_{a}x_{b}.
Proof.
To compute the sum of every element in a row, we fix i and u, and sum over j and v:
Σ_{j} Σ_{v} (M_{ij})_{uv}  =  Σ_{j} Σ_{v} (a[j, P_{j}[i,v]] b[i, Q_{i}[j,u]])^{n}  (by definition) 
 =  Σ_{j} b[i, Q_{i}[j,u]]^{n} Σ_{v} a[j, P_{j}[i,v]]^{n}  (since the b[] term doesn't depend on v) 
 =  Σ_{j} b[i, Q_{i}[j,u]]^{n} Σ_{v} a[j,v]^{n}  (since P_{j}[i,v] is just a permutation of v) 
 =  x_{a} Σ_{j} b[i, Q_{i}[j,u]]^{n}  (by the taxicab property for a) 
 =  x_{a} Σ_{j} b[i,j]^{n}  (since Q_{i}[j,u] is just a permutation of j) 
 =  x_{a}x_{b}  (by the taxicab property for b). 
So every row of the square M sums to x_{a}x_{b}. The proof that every column sums to x_{a}x_{b} is similar. End of proof.
Note that the sr elements in block M_{ij} 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 OnLine 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 = 1^{2} + 7^{2} = 5^{2} + 5^{2} can't be used to create a semimagic 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 semimagic square entries.
A 12×12 semimagic square of 4th powers
Let n = 4, r = 4, and s = 3.
I choose the taxicab(4,4,3)
650658  =  2^{4} + 16^{4} + 21^{4} + 25^{4} 
 =  5^{4} + 11^{4} + 12^{4} + 28^{4} 
 =  13^{4} + 19^{4} + 20^{4} + 24^{4} 
and the taxicab(4,3,4)
53809938  =  2^{4} + 71^{4} + 73^{4} 
 =  17^{4} + 62^{4} + 79^{4} 
 =  29^{4} + 53^{4} + 82^{4} 
 =  37^{4} + 46^{4} + 83^{4}. 
One can verify with a computer that all 144 products a_{ij}b_{uv} are distinct.
Let P_{1} = P_{2} = P_{3} =
Let Q_{1} = Q_{2} = Q_{3} = Q_{4} =
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 semimagic square
4^{4}  32^{4}  42^{4}  50^{4}  355^{4}  781^{4}  852^{4}  1988^{4}  949^{4}  1387^{4}  1460^{4}  1752^{4} 
142^{4}  1136^{4}  1491^{4}  1775^{4}  365^{4}  803^{4}  876^{4}  2044^{4}  26^{4}  38^{4}  40^{4}  48^{4} 
146^{4}  1168^{4}  1533^{4}  1825^{4}  10^{4}  22^{4}  24^{4}  56^{4}  923^{4}  1349^{4}  1420^{4}  1704^{4} 
272^{4}  357^{4}  425^{4}  34^{4}  682^{4}  744^{4}  1736^{4}  310^{4}  1501^{4}  1580^{4}  1896^{4}  1027^{4} 
992^{4}  1302^{4}  1550^{4}  124^{4}  869^{4}  948^{4}  2212^{4}  395^{4}  323^{4}  340^{4}  408^{4}  221^{4} 
1264^{4}  1659^{4}  1975^{4}  158^{4}  187^{4}  204^{4}  476^{4}  85^{4}  1178^{4}  1240^{4}  1488^{4}  806^{4} 
609^{4}  725^{4}  58^{4}  464^{4}  636^{4}  1484^{4}  265^{4}  583^{4}  1640^{4}  1968^{4}  1066^{4}  1558^{4} 
1113^{4}  1325^{4}  106^{4}  848^{4}  984^{4}  2296^{4}  410^{4}  902^{4}  580^{4}  696^{4}  377^{4}  551^{4} 
1722^{4}  2050^{4}  164^{4}  1312^{4}  348^{4}  812^{4}  145^{4}  319^{4}  1060^{4}  1272^{4}  689^{4}  1007^{4} 
925^{4}  74^{4}  592^{4}  777^{4}  1288^{4}  230^{4}  506^{4}  552^{4}  1992^{4}  1079^{4}  1577^{4}  1660^{4} 
1150^{4}  92^{4}  736^{4}  966^{4}  2324^{4}  415^{4}  913^{4}  996^{4}  888^{4}  481^{4}  703^{4}  740^{4} 
2075^{4}  166^{4}  1328^{4}  1743^{4}  1036^{4}  185^{4}  407^{4}  444^{4}  1104^{4}  598^{4}  874^{4}  920^{4} 
The magic sum is 650658 · 53809938 = 35011866639204. The largest entry is 28^{4}·83^{4} = 2324^{4}. 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 semimagic squares from other people
After I published my taxicab method formula on this webpage, Christian Boyer created:
Obtaining magic diagonals given a semimagic square
A random model
Semimagic squares produced by the taxicab method are unlikely to be fully magic. For a k × k semimagic 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^{2}, and c ≈ 1.867 for the 12×12 semimagic 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
p_{magic} = k/(2πc^{2}m^{2}),
or approximately 1/m^{2}. In the 12×12 semimagic 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 semimagic 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 5  0  1  2  3  4 
correction for k=1 
1.25  3.75  0.00  0.00  0.00 
correction for k=2 
0.31  1.88  2.81  0.00  0.00 
correction for k=3 
0.08  0.70  2.11  2.11  0.00 
... 
...  ...  ...  ...  ... 
correction for k=12 
1.22  0.83  0.67  0.97  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
f_{mod} = Π_{p prime} (correction factor when expressing m as a sum of k entries modulo p^{e})^{2}.
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 semimagic square above, f_{mod} ≈ 25.2, primarily due to effects modulo 5 and modulo 16.
Permuting rows and columns
Given a semimagic square, we can generate more semimagic 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
n_{perm} = (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×10^{12}.
 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+1x 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
p_{magic} · f_{mod} · n_{perm} =
f_{mod} k(k!)^{2} / (2πc^{2}m^{2} [k/2]! 2^{[k/2]+1})
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^{2})time algorithm
If n_{perm} > m^{2}, then a solution is likely to exist. We can try random row and column permutations. After ~m^{2} 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:
 Apply random permutations of rows (or columns) until we get a correct main diagonal. This should take ~m attempts.
 Generate a random permutation and apply it to both rows and columns until we get a correct antidiagonal. 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 semimagic square of 4th powers in October 2007. In his case, m=123784061778806 and k=44, so 44!/(22! 2^{22}) ≈ 5.64×10^{26} 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 antidiagonal is tricky as each set S and T should be symmetric (x is in S iff k+1x is in S). One can choose S = {1,2,...,q} union {k+1q, ..., k1, 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)!^{2}2^{k/2}) > m.
Obtaining magic diagonals given compatible taxicabs
In the previous section I discussed turning a given semimagic square into a fully magic square by permuting rows and columns. But if the semimagic square comes from the taxicab method, there is even more flexibility.
Rearranging the taxicabs
Note that in the taxicab(n,r,s)
x_{a}  =  a_{11}^{n} + a_{12}^{n} + ... + a_{1r}^{n} 
 =  a_{21}^{n} + a_{22}^{n} + ... + a_{2r}^{n} 
 =  ... 
 =  a_{s1}^{n} + a_{s2}^{n} + ... + a_{sr}^{n} 
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_{1}, ..., P_{s}, Q_{1}, ..., 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.
Order of the Latin square  Number of Latin squares with sorted first row 
1  1 
2  1 
3  2 
4  24 
5  1344 
6  1128960 
This gives us
n_{latin} = A000479(r)^{s}A000479(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
p_{magic} · f_{mod} · n_{perm} · n_{latin} = f_{mod} A000479(r)^{s}A000479(s)^{r} rs((rs)!)^{2} / (2πc^{2}m^{2} [rs/2]! 2^{[rs/2]+1}).
Search algorithms
The basic algorithm is:
 For every combination of Latin squares:
 Search the corresponding semimagic 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:

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  =  2^{4} + 21^{4} + 29^{4} + 32^{4} 
 =  7^{4} + 23^{4} + 24^{4} + 34^{4} 
 =  8^{4} + 9^{4} + 16^{4} + 37^{4} 
 =  14^{4} + 26^{4} + 27^{4} + 31^{4} 
and
321793923  =  3^{4} + 85^{4} + 97^{4} + 116^{4} 
 =  23^{4} + 25^{4} + 98^{4} + 123^{4} 
 =  43^{4} + 81^{4} + 95^{4} + 118^{4} 
 =  45^{4} + 73^{4} + 106^{4} + 113^{4}. 
One can verify with a computer that all 256 products a_{ij}b_{uv} are distinct. As we know by now, this gives a 16×16 semimagic square with magic sum m = 1950354 · 321793923 = 627612064898742. The largest entry is 37^{4}·123^{4} = 4551^{4}. 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^{8} 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
 their sum equals the magic sum, and
 their positions can potentially form a diagonal after rearrangement.
To explain what I mean by (2), recall from the taxicab method (semimagic 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 14  cols 58  cols 912  cols 1316  sum 
rows 14  0  1  2  1  4 
rows 58  2  1  1  0  4 
rows 912  2  1  0  1  4 
rows 1316  0  1  1  2  4 
sum  4  4  4  4  16 
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(p_{magic} · f_{mod}) ≈ 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.
6^{4}  63^{4}  87^{4}  96^{4} 
170^{4}  1785^{4}  2465^{4}  2720^{4} 
194^{4}  2037^{4}  2813^{4}  3104^{4} 
232^{4}  2436^{4}  3364^{4}  3712^{4} 

679^{4}  2231^{4}  2328^{4}  3298^{4} 
21^{4}  69^{4}  72^{4}  102^{4} 
812^{4}  2668^{4}  2784^{4}  3944^{4} 
595^{4}  1955^{4}  2040^{4}  2890^{4} 

928^{4}  1044^{4}  1856^{4}  4292^{4} 
776^{4}  873^{4}  1552^{4}  3589^{4} 
680^{4}  765^{4}  1360^{4}  3145^{4} 
24^{4}  27^{4}  48^{4}  111^{4} 

1190^{4}  2210^{4}  2295^{4}  2635^{4} 
1624^{4}  3016^{4}  3132^{4}  3596^{4} 
42^{4}  78^{4}  81^{4}  93^{4} 
1358^{4}  2522^{4}  2619^{4}  3007^{4} 

736^{4}  667^{4}  483^{4}  46^{4} 
800^{4}  725^{4}  525^{4}  50^{4} 
3136^{4}  2842^{4}  2058^{4}  196^{4} 
3936^{4}  3567^{4}  2583^{4}  246^{4} 

4182^{4}  861^{4}  2829^{4}  2952^{4} 
3332^{4}  686^{4}  2254^{4}  2352^{4} 
782^{4}  161^{4}  529^{4}  552^{4} 
850^{4}  175^{4}  575^{4}  600^{4} 

400^{4}  200^{4}  925^{4}  225^{4} 
368^{4}  184^{4}  851^{4}  207^{4} 
1968^{4}  984^{4}  4551^{4}  1107^{4} 
1568^{4}  784^{4}  3626^{4}  882^{4} 

3038^{4}  2646^{4}  1372^{4}  2548^{4} 
3813^{4}  3321^{4}  1722^{4}  3198^{4} 
775^{4}  675^{4}  350^{4}  650^{4} 
713^{4}  621^{4}  322^{4}  598^{4} 

903^{4}  86^{4}  1376^{4}  1247^{4} 
1701^{4}  162^{4}  2592^{4}  2349^{4} 
1995^{4}  190^{4}  3040^{4}  2755^{4} 
2478^{4}  236^{4}  3776^{4}  3422^{4} 

2185^{4}  2280^{4}  3230^{4}  665^{4} 
989^{4}  1032^{4}  1462^{4}  301^{4} 
2714^{4}  2832^{4}  4012^{4}  826^{4} 
1863^{4}  1944^{4}  2754^{4}  567^{4} 

729^{4}  2997^{4}  648^{4}  1296^{4} 
1062^{4}  4366^{4}  944^{4}  1888^{4} 
387^{4}  1591^{4}  344^{4}  688^{4} 
855^{4}  3515^{4}  760^{4}  1520^{4} 

3068^{4}  1652^{4}  3658^{4}  3186^{4} 
2470^{4}  1330^{4}  2945^{4}  2565^{4} 
2106^{4}  1134^{4}  2511^{4}  2187^{4} 
1118^{4}  602^{4}  1333^{4}  1161^{4} 

1305^{4}  1440^{4}  90^{4}  945^{4} 
2117^{4}  2336^{4}  146^{4}  1533^{4} 
3074^{4}  3392^{4}  212^{4}  2226^{4} 
3277^{4}  3616^{4}  226^{4}  2373^{4} 

1752^{4}  2482^{4}  511^{4}  1679^{4} 
1080^{4}  1530^{4}  315^{4}  1035^{4} 
2712^{4}  3842^{4}  791^{4}  2599^{4} 
2544^{4}  3604^{4}  742^{4}  2438^{4} 

4181^{4}  1808^{4}  1017^{4}  904^{4} 
3922^{4}  1696^{4}  954^{4}  848^{4} 
2701^{4}  1168^{4}  657^{4}  584^{4} 
1665^{4}  720^{4}  405^{4}  360^{4} 

2862^{4}  3286^{4}  2756^{4}  1484^{4} 
3051^{4}  3503^{4}  2938^{4}  1582^{4} 
1215^{4}  1395^{4}  1170^{4}  630^{4} 
1971^{4}  2263^{4}  1898^{4}  1022^{4} 

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.
3298^{4}  928^{4}  1044^{4}  2295^{4}  63^{4}  96^{4}  2231^{4}  679^{4}  2210^{4}  87^{4}  4292^{4}  2635^{4}  6^{4}  2328^{4}  1856^{4}  1190^{4} 
102^{4}  776^{4}  873^{4}  3132^{4}  1785^{4}  2720^{4}  69^{4}  21^{4}  3016^{4}  2465^{4}  3589^{4}  3596^{4}  170^{4}  72^{4}  1552^{4}  1624^{4} 
3944^{4}  680^{4}  765^{4}  81^{4}  2037^{4}  3104^{4}  2668^{4}  812^{4}  78^{4}  2813^{4}  3145^{4}  93^{4}  194^{4}  2784^{4}  1360^{4}  42^{4} 
2890^{4}  24^{4}  27^{4}  2619^{4}  2436^{4}  3712^{4}  1955^{4}  595^{4}  2522^{4}  3364^{4}  111^{4}  3007^{4}  232^{4}  2040^{4}  48^{4}  1358^{4} 
1035^{4}  3922^{4}  1696^{4}  2938^{4}  2336^{4}  1533^{4}  1530^{4}  1080^{4}  3503^{4}  146^{4}  848^{4}  1582^{4}  2117^{4}  315^{4}  954^{4}  3051^{4} 
2352^{4}  368^{4}  184^{4}  1722^{4}  725^{4}  50^{4}  686^{4}  3332^{4}  3321^{4}  525^{4}  207^{4}  3198^{4}  800^{4}  2254^{4}  851^{4}  3813^{4} 
552^{4}  1968^{4}  984^{4}  350^{4}  2842^{4}  196^{4}  161^{4}  782^{4}  675^{4}  2058^{4}  1107^{4}  650^{4}  3136^{4}  529^{4}  4551^{4}  775^{4} 
2599^{4}  2701^{4}  1168^{4}  1170^{4}  3392^{4}  2226^{4}  3842^{4}  2712^{4}  1395^{4}  212^{4}  584^{4}  630^{4}  3074^{4}  791^{4}  657^{4}  1215^{4} 
301^{4}  1062^{4}  4366^{4}  2945^{4}  162^{4}  2349^{4}  1032^{4}  989^{4}  1330^{4}  2592^{4}  1888^{4}  2565^{4}  1701^{4}  1462^{4}  944^{4}  2470^{4} 
567^{4}  855^{4}  3515^{4}  1333^{4}  236^{4}  3422^{4}  1944^{4}  1863^{4}  602^{4}  3776^{4}  1520^{4}  1161^{4}  2478^{4}  2754^{4}  760^{4}  1118^{4} 
600^{4}  1568^{4}  784^{4}  322^{4}  3567^{4}  246^{4}  175^{4}  850^{4}  621^{4}  2583^{4}  882^{4}  598^{4}  3936^{4}  575^{4}  3626^{4}  713^{4} 
2438^{4}  1665^{4}  720^{4}  1898^{4}  3616^{4}  2373^{4}  3604^{4}  2544^{4}  2263^{4}  226^{4}  360^{4}  1022^{4}  3277^{4}  742^{4}  405^{4}  1971^{4} 
826^{4}  387^{4}  1591^{4}  2511^{4}  190^{4}  2755^{4}  2832^{4}  2714^{4}  1134^{4}  3040^{4}  688^{4}  2187^{4}  1995^{4}  4012^{4}  344^{4}  2106^{4} 
2952^{4}  400^{4}  200^{4}  1372^{4}  667^{4}  46^{4}  861^{4}  4182^{4}  2646^{4}  483^{4}  225^{4}  2548^{4}  736^{4}  2829^{4}  925^{4}  3038^{4} 
1679^{4}  4181^{4}  1808^{4}  2756^{4}  1440^{4}  945^{4}  2482^{4}  1752^{4}  3286^{4}  90^{4}  904^{4}  1484^{4}  1305^{4}  511^{4}  1017^{4}  2862^{4} 
665^{4}  729^{4}  2997^{4}  3658^{4}  86^{4}  1247^{4}  2280^{4}  2185^{4}  1652^{4}  1376^{4}  1296^{4}  3186^{4}  903^{4}  3230^{4}  648^{4}  3068^{4} 
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 semimagic square of 4th powers, I applied the same method above to find magic diagonals.
Here's Christian's published semimagic square, with the elements that will become diagonals highlighted:
2^{4}  10^{4}  24^{4}  32^{4}  58^{4} 
109^{4}  545^{4}  1308^{4}  1744^{4}  3161^{4} 
111^{4}  555^{4}  1332^{4}  1776^{4}  3219^{4} 

763^{4}  1090^{4}  2071^{4}  2289^{4}  2834^{4} 
777^{4}  1110^{4}  2109^{4}  2331^{4}  2886^{4} 
14^{4}  20^{4}  38^{4}  42^{4}  52^{4} 

888^{4}  1221^{4}  1998^{4}  2553^{4}  2775^{4} 
16^{4}  22^{4}  36^{4}  46^{4}  50^{4} 
872^{4}  1199^{4}  1962^{4}  2507^{4}  2725^{4} 

105^{4}  252^{4}  336^{4}  609^{4}  21^{4} 
490^{4}  1176^{4}  1568^{4}  2842^{4}  98^{4} 
595^{4}  1428^{4}  1904^{4}  3451^{4}  119^{4} 

980^{4}  1862^{4}  2058^{4}  2548^{4}  686^{4} 
1190^{4}  2261^{4}  2499^{4}  3094^{4}  833^{4} 
210^{4}  399^{4}  441^{4}  546^{4}  147^{4} 

1309^{4}  2142^{4}  2737^{4}  2975^{4}  952^{4} 
231^{4}  378^{4}  483^{4}  525^{4}  168^{4} 
1078^{4}  1764^{4}  2254^{4}  2450^{4}  784^{4} 

324^{4}  432^{4}  783^{4}  27^{4}  135^{4} 
1128^{4}  1504^{4}  2726^{4}  94^{4}  470^{4} 
1452^{4}  1936^{4}  3509^{4}  121^{4}  605^{4} 

1786^{4}  1974^{4}  2444^{4}  658^{4}  940^{4} 
2299^{4}  2541^{4}  3146^{4}  847^{4}  1210^{4} 
513^{4}  567^{4}  702^{4}  189^{4}  270^{4} 

2178^{4}  2783^{4}  3025^{4}  968^{4}  1331^{4} 
486^{4}  621^{4}  675^{4}  216^{4}  297^{4} 
1692^{4}  2162^{4}  2350^{4}  752^{4}  1034^{4} 

544^{4}  986^{4}  34^{4}  170^{4}  408^{4} 
1424^{4}  2581^{4}  89^{4}  445^{4}  1068^{4} 
1968^{4}  3567^{4}  123^{4}  615^{4}  1476^{4} 

1869^{4}  2314^{4}  623^{4}  890^{4}  1691^{4} 
2583^{4}  3198^{4}  861^{4}  1230^{4}  2337^{4} 
714^{4}  884^{4}  238^{4}  340^{4}  646^{4} 

2829^{4}  3075^{4}  984^{4}  1353^{4}  2214^{4} 
782^{4}  850^{4}  272^{4}  374^{4}  612^{4} 
2047^{4}  2225^{4}  712^{4}  979^{4}  1602^{4} 

1769^{4}  61^{4}  305^{4}  732^{4}  976^{4} 
1914^{4}  66^{4}  330^{4}  792^{4}  1056^{4} 
3683^{4}  127^{4}  635^{4}  1524^{4}  2032^{4} 

1716^{4}  462^{4}  660^{4}  1254^{4}  1386^{4} 
3302^{4}  889^{4}  1270^{4}  2413^{4}  2667^{4} 
1586^{4}  427^{4}  610^{4}  1159^{4}  1281^{4} 

3175^{4}  1016^{4}  1397^{4}  2286^{4}  2921^{4} 
1525^{4}  488^{4}  671^{4}  1098^{4}  1403^{4} 
1650^{4}  528^{4}  726^{4}  1188^{4}  1518^{4} 

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.
2^{4}  10^{4}  24^{4}  32^{4}  58^{4} 
109^{4}  545^{4}  1308^{4}  1744^{4}  3161^{4} 
111^{4}  555^{4}  1332^{4}  1776^{4}  3219^{4} 

763^{4}  1090^{4}  2071^{4}  2289^{4}  2834^{4} 
777^{4}  1110^{4}  2109^{4}  2331^{4}  2886^{4} 
14^{4}  20^{4}  38^{4}  42^{4}  52^{4} 

888^{4}  1221^{4}  1998^{4}  2553^{4}  2775^{4} 
16^{4}  22^{4}  36^{4}  46^{4}  50^{4} 
872^{4}  1199^{4}  1962^{4}  2507^{4}  2725^{4} 

105^{4}  252^{4}  336^{4}  609^{4}  21^{4} 
490^{4}  1176^{4}  1568^{4}  2842^{4}  98^{4} 
595^{4}  1428^{4}  1904^{4}  3451^{4}  119^{4} 

980^{4}  686^{4}  2058^{4}  2548^{4}  1862^{4} 
1190^{4}  833^{4}  2499^{4}  3094^{4}  2261^{4} 
210^{4}  147^{4}  441^{4}  546^{4}  399^{4} 

1309^{4}  2737^{4}  2975^{4}  2142^{4}  952^{4} 
231^{4}  483^{4}  525^{4}  378^{4}  168^{4} 
1078^{4}  2254^{4}  2450^{4}  1764^{4}  784^{4} 

783^{4}  432^{4}  27^{4}  324^{4}  135^{4} 
2726^{4}  1504^{4}  94^{4}  1128^{4}  470^{4} 
3509^{4}  1936^{4}  121^{4}  1452^{4}  605^{4} 

2541^{4}  2299^{4}  3146^{4}  1210^{4}  847^{4} 
567^{4}  513^{4}  702^{4}  270^{4}  189^{4} 
1974^{4}  1786^{4}  2444^{4}  940^{4}  658^{4} 

1692^{4}  2350^{4}  1034^{4}  752^{4}  2162^{4} 
2178^{4}  3025^{4}  1331^{4}  968^{4}  2783^{4} 
486^{4}  675^{4}  297^{4}  216^{4}  621^{4} 

544^{4}  34^{4}  986^{4}  170^{4}  408^{4} 
1424^{4}  89^{4}  2581^{4}  445^{4}  1068^{4} 
1968^{4}  123^{4}  3567^{4}  615^{4}  1476^{4} 

1691^{4}  2314^{4}  890^{4}  623^{4}  1869^{4} 
2337^{4}  3198^{4}  1230^{4}  861^{4}  2583^{4} 
646^{4}  884^{4}  340^{4}  238^{4}  714^{4} 

3075^{4}  984^{4}  2829^{4}  1353^{4}  2214^{4} 
850^{4}  272^{4}  782^{4}  374^{4}  612^{4} 
2225^{4}  712^{4}  2047^{4}  979^{4}  1602^{4} 

732^{4}  1769^{4}  305^{4}  61^{4}  976^{4} 
792^{4}  1914^{4}  330^{4}  66^{4}  1056^{4} 
1524^{4}  3683^{4}  635^{4}  127^{4}  2032^{4} 

1716^{4}  1386^{4}  462^{4}  1254^{4}  660^{4} 
3302^{4}  2667^{4}  889^{4}  2413^{4}  1270^{4} 
1586^{4}  1281^{4}  427^{4}  1159^{4}  610^{4} 

2921^{4}  2286^{4}  1016^{4}  3175^{4}  1397^{4} 
1403^{4}  1098^{4}  488^{4}  1525^{4}  671^{4} 
1518^{4}  1188^{4}  528^{4}  1650^{4}  726^{4} 

Here's the fully magic square, after permuting the rows and columns:
1424^{4}  89^{4}  2581^{4}  374^{4}  272^{4}  2337^{4}  3198^{4}  850^{4}  782^{4}  1230^{4}  612^{4}  1068^{4}  2583^{4}  445^{4}  861^{4} 
3509^{4}  1936^{4}  121^{4}  216^{4}  675^{4}  1974^{4}  1786^{4}  486^{4}  297^{4}  2444^{4}  621^{4}  605^{4}  658^{4}  1452^{4}  940^{4} 
109^{4}  545^{4}  1308^{4}  46^{4}  22^{4}  777^{4}  1110^{4}  16^{4}  36^{4}  2109^{4}  50^{4}  3161^{4}  2886^{4}  1744^{4}  2331^{4} 
544^{4}  34^{4}  986^{4}  1353^{4}  984^{4}  1691^{4}  2314^{4}  3075^{4}  2829^{4}  890^{4}  2214^{4}  408^{4}  1869^{4}  170^{4}  623^{4} 
2^{4}  10^{4}  24^{4}  2553^{4}  1221^{4}  763^{4}  1090^{4}  888^{4}  1998^{4}  2071^{4}  2775^{4}  58^{4}  2834^{4}  32^{4}  2289^{4} 
732^{4}  1769^{4}  305^{4}  3175^{4}  2286^{4}  1716^{4}  1386^{4}  2921^{4}  1016^{4}  462^{4}  1397^{4}  976^{4}  660^{4}  61^{4}  1254^{4} 
792^{4}  1914^{4}  330^{4}  1525^{4}  1098^{4}  3302^{4}  2667^{4}  1403^{4}  488^{4}  889^{4}  671^{4}  1056^{4}  1270^{4}  66^{4}  2413^{4} 
1968^{4}  123^{4}  3567^{4}  979^{4}  712^{4}  646^{4}  884^{4}  2225^{4}  2047^{4}  340^{4}  1602^{4}  1476^{4}  714^{4}  615^{4}  238^{4} 
490^{4}  1176^{4}  1568^{4}  378^{4}  483^{4}  1190^{4}  833^{4}  231^{4}  525^{4}  2499^{4}  168^{4}  98^{4}  2261^{4}  2842^{4}  3094^{4} 
783^{4}  432^{4}  27^{4}  752^{4}  2350^{4}  2541^{4}  2299^{4}  1692^{4}  1034^{4}  3146^{4}  2162^{4}  135^{4}  847^{4}  324^{4}  1210^{4} 
595^{4}  1428^{4}  1904^{4}  1764^{4}  2254^{4}  210^{4}  147^{4}  1078^{4}  2450^{4}  441^{4}  784^{4}  119^{4}  399^{4}  3451^{4}  546^{4} 
2726^{4}  1504^{4}  94^{4}  968^{4}  3025^{4}  567^{4}  513^{4}  2178^{4}  1331^{4}  702^{4}  2783^{4}  470^{4}  189^{4}  1128^{4}  270^{4} 
105^{4}  252^{4}  336^{4}  2142^{4}  2737^{4}  980^{4}  686^{4}  1309^{4}  2975^{4}  2058^{4}  952^{4}  21^{4}  1862^{4}  609^{4}  2548^{4} 
111^{4}  555^{4}  1332^{4}  2507^{4}  1199^{4}  14^{4}  20^{4}  872^{4}  1962^{4}  38^{4}  2725^{4}  3219^{4}  52^{4}  1776^{4}  42^{4} 
1524^{4}  3683^{4}  635^{4}  1650^{4}  1188^{4}  1586^{4}  1281^{4}  1518^{4}  528^{4}  427^{4}  726^{4}  2032^{4}  610^{4}  127^{4}  1159^{4} 
The magic sum is 794179 · 292965218 = 232666823866022. The largest entry is 29^{4}·127^{4} = 3683^{4}. 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 semimagic squares that could lead to new records if their entries can be rearranged to form 2 magic diagonals:
 Christian Boyer's 25x25 semimagic square of 5th powers,
 the 25×25 semimagic square of 5th powers with the smallest possible magic sum, using the two taxicabs(5,5,5)
1189260043200  =  22^{5} + 58^{5} + 74^{5} + 161^{5} + 255^{5} 
 =  28^{5} + 62^{5} + 155^{5} + 165^{5} + 250^{5} 
 =  43^{5} + 149^{5} + 172^{5} + 176^{5} + 240^{5} 
 =  79^{5} + 119^{5} + 186^{5} + 195^{5} + 231^{5} 
 =  123^{5} + 134^{5} + 187^{5} + 205^{5} + 221^{5} 
and
39931433836700  =  10^{5} + 223^{5} + 292^{5} + 382^{5} + 493^{5} 
 =  33^{5} + 191^{5} + 299^{5} + 385^{5} + 492^{5} 
 =  35^{5} + 199^{5} + 271^{5} + 412^{5} + 483^{5} 
 =  113^{5} + 227^{5} + 365^{5} + 416^{5} + 459^{5} 
 =  263^{5} + 332^{5} + 389^{5} + 406^{5} + 430^{5}. 
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 semimagic squares I discussed earlier, sorted by power and by size.
power  semimagic square  magic sum  c  f_{mod}  p_{magic} · f_{mod} · n_{perm}  p_{magic} · f_{mod} · n_{perm} · n_{latin}  magic square 
4  12×12 FL  3.50×10^{13}  1.87  25.2  2.8×10^{14}  6.2×10^{9}  none (FL) 
15×15 CB  2.33×10^{14}  2.00  17.7  2.6×10^{10}  20  found (FL) 
16×16 FL  6.28×10^{14}  1.78  14.1  6.1×10^{10}  67  found (FL) 
5  25×25 FL  4.75×10^{25}  2.32  1.22  2.5×10^{14}  4.7×10^{17}  ? 
25×25 CB  8.28×10^{25}  2.65  3.02  1.5×10^{14}  2.9×10^{17}  ? 
The results show that a 25×25 magic square of 5th powers can almost certainly be obtained from the respective semimagic 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 semimagic square of nth powers by using three taxicab(n,2,2) numbers. This hints at a general method to construct a (rst)×(rst) semimagic 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 semimagic 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