There are only seven (six unique) solutions up to row/column sums of $10^8$ (of course there are 1152 equivalent squares by permutation/transposition for each one):
$$
\left(
\begin{array}{cccc}
2 & 16 & 108 & 180 \\
24 & 192 & 9 & 15 \\
144 & 18 & 150 & 90 \\
160 & 20 & 135 & 81 \\
\end{array}
\right),\left(
\begin{array}{cccc}
10 & 27 & 228 & 288 \\
120 & 324 & 19 & 24 \\
243 & 90 & 240 & 190 \\
270 & 100 & 216 & 171 \\
\end{array}
\right),\left(
\begin{array}{cccc}
25 & 207 & 243 & 269 \\
239 & 135 & 297 & 73 \\
271 & 89 & 109 & 275 \\
209 & 313 & 95 & 127 \\
\end{array}
\right),\\
\left(
\begin{array}{cccc}
4 & 32 & 216 & 360 \\
48 & 384 & 18 & 30 \\
288 & 36 & 300 & 180 \\
320 & 40 & 270 & 162 \\
\end{array}
\right),\left(
\begin{array}{cccc}
2 & 24 & 306 & 340 \\
180 & 15 & 330 & 297 \\
34 & 408 & 18 & 20 \\
396 & 33 & 150 & 135 \\
\end{array}
\right),\left(
\begin{array}{cccc}
9 & 34 & 192 & 396 \\
108 & 408 & 16 & 33 \\
306 & 81 & 330 & 160 \\
340 & 90 & 297 & 144 \\
\end{array}
\right),\\
\left(
\begin{array}{cccc}
20 & 54 & 304 & 384 \\
160 & 432 & 38 & 48 \\
243 & 90 & 360 & 285 \\
405 & 150 & 216 & 171 \\
\end{array}
\right)
$$
The fourth solution is twice the first one.
I wrote a Mathematica code that searches up to a specific sum, but was not able to search beyond $300\,000$; only after translating the code to C I could find the above solutions (the smallest one having a row/column sum of $7\,095\,816$). Any Mathematica code will have to be very efficient to search beyond sums of 7 million, and I think Mathematica is not the right tool for this job.
For instance, the first solution has 65 decompositions into a sum of four cubes:
PowersRepresentations[7095816, 4, 3]
(* {{1, 23, 76, 188}, {1, 62, 90, 183}, {2, 10, 118, 176}, {2, 16, 24, 192},
{2, 16, 108, 180}, {2, 24, 144, 160}, {4, 34, 146, 158}, {4, 50, 110, 178},
{6, 11, 70, 189}, {6, 11, 133, 168}, {6, 60, 138, 162}, {6, 64, 89, 183},
{6, 66, 72, 186}, {6, 114, 138, 144}, {8, 30, 145, 159}, {8, 68, 82, 184},
{9, 12, 87, 186}, {9, 15, 24, 192}, {9, 15, 108, 180}, {9, 108, 135, 150},
{10, 23, 113, 178}, {12, 25, 48, 191}, {14, 36, 147, 157}, {14, 46, 77, 187},
{15, 45, 63, 189}, {15, 66, 108, 177}, {15, 81, 90, 180}, {16, 18, 20, 192},
{16, 101, 102, 171}, {18, 20, 144, 160}, {18, 72, 96, 180}, {18, 90, 144, 150},
{20, 36, 97, 183}, {20, 81, 135, 160}, {24, 80, 114, 172}, {24, 90, 123, 165},
{27, 44, 102, 181}, {27, 109, 134, 150}, {31, 86, 96, 177}, {33, 36, 96, 183},
{34, 57, 135, 164}, {36, 52, 146, 156}, {36, 82, 144, 152}, {37, 58, 59, 188},
{40, 50, 64, 188}, {42, 84, 121, 167}, {43, 48, 120, 173}, {44, 48, 99, 181},
{49, 84, 120, 167}, {50, 67, 101, 178}, {50, 82, 110, 172}, {54, 60, 66, 186},
{54, 86, 120, 166}, {54, 124, 135, 137}, {62, 66, 114, 172}, {62, 79, 97, 176},
{64, 104, 130, 152}, {71, 102, 138, 145}, {72, 90, 139, 149}, {75, 96, 141, 144},
{76, 81, 95, 174}, {81, 90, 135, 150}, {84, 91, 125, 156}, {86, 90, 123, 157},
{94, 102, 108, 158}} *)
C code
Ugly but fast code that finds all tuples $\{i_0,i_1,i_2,i_3\}$ whose cubes sum to a given value, and then tries all permutations to see if we can get a cube-magic square. The checks and permutations are hard-coded (using Mathematica to generate some code, at least).
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <inttypes.h>
#define MAXSUMS 100
typedef struct {
int number_of_sums;
uint32_t sum[4*MAXSUMS];
} four_sums;
static inline void
check(const uint32_t a00, const uint32_t a01, const uint32_t a02, const uint32_t a03,
const uint32_t a10, const uint32_t a11, const uint32_t a12, const uint32_t a13,
const uint32_t a20, const uint32_t a21, const uint32_t a22, const uint32_t a23,
const uint32_t a30, const uint32_t a31, const uint32_t a32, const uint32_t a33)
{
uint32_t s0 = a00+a10+a20+a30;
uint32_t s1 = a01+a11+a21+a31;
uint32_t s2 = a02+a12+a22+a32;
uint32_t s3 = a03+a13+a23+a33;
if (s0 == s1 && s0 == s2 && s0 == s3) {
printf("{{%"PRIu32",%"PRIu32",%"PRIu32",%"PRIu32"},"
"{%"PRIu32",%"PRIu32",%"PRIu32",%"PRIu32"},"
"{%"PRIu32",%"PRIu32",%"PRIu32",%"PRIu32"},"
"{%"PRIu32",%"PRIu32",%"PRIu32",%"PRIu32"}}\n",
a00,a01,a02,a03,a10,a11,a12,a13,a20,a21,a22,a23,a30,a31,a32,a33);
}
}
static inline void
permute3(const uint32_t a00, const uint32_t a01, const uint32_t a02, const uint32_t a03,
const uint32_t a10, const uint32_t a11, const uint32_t a12, const uint32_t a13,
const uint32_t a20, const uint32_t a21, const uint32_t a22, const uint32_t a23,
const uint32_t a30, const uint32_t a31, const uint32_t a32, const uint32_t a33)
{
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a31,a32,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a31,a33,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a32,a31,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a32,a33,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a33,a31,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a33,a32,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a30,a32,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a30,a33,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a32,a30,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a32,a33,a30);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a33,a30,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a31,a33,a32,a30);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a30,a31,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a30,a33,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a31,a30,a33);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a31,a33,a30);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a33,a30,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a32,a33,a31,a30);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a30,a31,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a30,a32,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a31,a30,a32);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a31,a32,a30);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a32,a30,a31);
check(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a33,a32,a31,a30);
}
static inline void
permute2(const uint32_t a00, const uint32_t a01, const uint32_t a02, const uint32_t a03,
const uint32_t a10, const uint32_t a11, const uint32_t a12, const uint32_t a13,
const uint32_t a20, const uint32_t a21, const uint32_t a22, const uint32_t a23,
const uint32_t a30, const uint32_t a31, const uint32_t a32, const uint32_t a33)
{
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a23,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a22,a21,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a22,a23,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a23,a21,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a20,a23,a22,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a20,a22,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a20,a23,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a22,a20,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a22,a23,a20, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a23,a20,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a21,a23,a22,a20, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a20,a21,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a20,a23,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a21,a20,a23, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a21,a23,a20, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a23,a20,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a22,a23,a21,a20, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a20,a21,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a20,a22,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a21,a20,a22, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a21,a22,a20, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a22,a20,a21, a30,a31,a32,a33);
permute3(a00,a01,a02,a03, a10,a11,a12,a13, a23,a22,a21,a20, a30,a31,a32,a33);
}
static inline void
permute1(const uint32_t a00, const uint32_t a01, const uint32_t a02, const uint32_t a03,
const uint32_t a10, const uint32_t a11, const uint32_t a12, const uint32_t a13,
const uint32_t a20, const uint32_t a21, const uint32_t a22, const uint32_t a23,
const uint32_t a30, const uint32_t a31, const uint32_t a32, const uint32_t a33)
{
permute2(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a10,a11,a13,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a10,a12,a11,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a10,a12,a13,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a10,a13,a11,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a10,a13,a12,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a10,a12,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a10,a13,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a12,a10,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a12,a13,a10, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a13,a10,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a11,a13,a12,a10, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a10,a11,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a10,a13,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a11,a10,a13, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a11,a13,a10, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a13,a10,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a12,a13,a11,a10, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a10,a11,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a10,a12,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a11,a10,a12, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a11,a12,a10, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a12,a10,a11, a20,a21,a22,a23, a30,a31,a32,a33);
permute2(a00,a01,a02,a03, a13,a12,a11,a10, a20,a21,a22,a23, a30,a31,a32,a33);
}
static void
check_sums(const uint32_t a00, const uint32_t a01, const uint32_t a02, const uint32_t a03,
const uint32_t a10, const uint32_t a11, const uint32_t a12, const uint32_t a13,
const uint32_t a20, const uint32_t a21, const uint32_t a22, const uint32_t a23,
const uint32_t a30, const uint32_t a31, const uint32_t a32, const uint32_t a33)
{
if (a00!=a10 && a00!=a11 && a00!=a12 && a00!=a13 && a00!=a20 && a00!=a21 &&
a00!=a22 && a00!=a23 && a00!=a30 && a00!=a31 && a00!=a32 && a00!=a33 &&
a01!=a10 && a01!=a11 && a01!=a12 && a01!=a13 && a01!=a20 && a01!=a21 &&
a01!=a22 && a01!=a23 && a01!=a30 && a01!=a31 && a01!=a32 && a01!=a33 &&
a02!=a10 && a02!=a11 && a02!=a12 && a02!=a13 && a02!=a20 && a02!=a21 &&
a02!=a22 && a02!=a23 && a02!=a30 && a02!=a31 && a02!=a32 && a02!=a33 &&
a03!=a10 && a03!=a11 && a03!=a12 && a03!=a13 && a03!=a20 && a03!=a21 &&
a03!=a22 && a03!=a23 && a03!=a30 && a03!=a31 && a03!=a32 && a03!=a33 &&
a10!=a20 && a10!=a21 && a10!=a22 && a10!=a23 && a10!=a30 && a10!=a31 &&
a10!=a32 && a10!=a33 && a11!=a20 && a11!=a21 && a11!=a22 && a11!=a23 &&
a11!=a30 && a11!=a31 && a11!=a32 && a11!=a33 && a12!=a20 && a12!=a21 &&
a12!=a22 && a12!=a23 && a12!=a30 && a12!=a31 && a12!=a32 && a12!=a33 &&
a13!=a20 && a13!=a21 && a13!=a22 && a13!=a23 && a13!=a30 && a13!=a31 &&
a13!=a32 && a13!=a33 && a20!=a30 && a20!=a31 && a20!=a32 && a20!=a33 &&
a21!=a30 && a21!=a31 && a21!=a32 && a21!=a33 && a22!=a30 && a22!=a31 &&
a22!=a32 && a22!=a33 && a23!=a30 && a23!=a31 && a23!=a32 && a23!=a33)
permute1(a00,a01,a02,a03, a10,a11,a12,a13, a20,a21,a22,a23, a30,a31,a32,a33);
}
int four_cubes(const uint32_t m)
{
// compute the cube-root of the upper limit m
uint32_t m3 = pow(m+1, 1./3);
// compute the cubes of all numbers from 0 to m3
uint32_t *cube = malloc((m3+1)*sizeof(uint32_t));
if (cube == NULL) {
fprintf(stderr, "%s: malloc failure\n", __func__);
return 3;
}
for (uint32_t i=0; i<=m3; i++)
cube[i] = i*i*i;
// compute all possible decompositions of the numbers 0..m into four cubes
four_sums *F = malloc((m+1)*sizeof(four_sums));
if (F == NULL) {
fprintf(stderr, "%s: malloc failure\n", __func__);
return 3;
}
for (uint32_t i=0; i<=m; i++)
F[i].number_of_sums = 0;
for (uint32_t i0=0; i0<=m3; i0++)
for (uint32_t i1=i0+1; i1<=m3; i1++)
for (uint32_t i2=i1+1; i2<=m3; i2++)
for (uint32_t i3=i2+1; i3<=m3; i3++) {
uint32_t s = cube[i0]+cube[i1]+cube[i2]+cube[i3];
if (s <= m) {
if (F[s].number_of_sums == MAXSUMS) {
fprintf(stderr, "exceeded MAXSUMS=%d at s=%"PRIu32"\n", MAXSUMS, s);
return 4;
}
F[s].sum[4*F[s].number_of_sums+0] = cube[i0];
F[s].sum[4*F[s].number_of_sums+1] = cube[i1];
F[s].sum[4*F[s].number_of_sums+2] = cube[i2];
F[s].sum[4*F[s].number_of_sums+3] = cube[i3];
F[s].number_of_sums++;
}
}
// check if any of the cube-decompositions gives a perfect square
for (uint32_t i=0; i<=m; i++)
if (F[i].number_of_sums >= 8)
for (int i0=0; i0<F[i].number_of_sums; i0++)
for (int i1=i0+1; i1<F[i].number_of_sums; i1++)
for (int i2=i1+1; i2<F[i].number_of_sums; i2++)
for (int i3=i2+1; i3<F[i].number_of_sums; i3++)
check_sums(F[i].sum[4*i0+0], F[i].sum[4*i0+1], F[i].sum[4*i0+2], F[i].sum[4*i0+3],
F[i].sum[4*i1+0], F[i].sum[4*i1+1], F[i].sum[4*i1+2], F[i].sum[4*i1+3],
F[i].sum[4*i2+0], F[i].sum[4*i2+1], F[i].sum[4*i2+2], F[i].sum[4*i2+3],
F[i].sum[4*i3+0], F[i].sum[4*i3+1], F[i].sum[4*i3+2], F[i].sum[4*i3+3]);
free(cube);
free(F);
return 0;
}
void show_usage(char const * const name)
{
fprintf(stderr, "usage: %s <max>\n", name);
fprintf(stderr, "Look for 4x4 squares of integers that have the same cube-sum"
" in every row and column.\n");
fprintf(stderr, "Scan numbers up to <max>.\n");
}
int main(int argc, char *argv[]) {
if (argc != 2) {
show_usage(argv[0]);
return 1;
}
long int m;
if (sscanf(argv[1], "%ld", &m) != 1) {
show_usage(argv[0]);
return 2;
}
return four_cubes(m);
}
Compile with
gcc -Ofast -std=c99 -lm fourcubes.c -o fourcubes
and run with
./fourcubes 10000000
Update
Here is a list of all 39 solutions with cubic row/column sums up to $10^9$. Maybe someone can find some order in these numbers.
{{{2,16,108,180},{24,192,9,15},{144,18,150,90},{160,20,135,81}},
{{10,27,228,288},{120,324,19,24},{243,90,240,190},{270,100,216,171}},
{{25,207,243,269},{239,135,297,73},{271,89,109,275},{209,313,95,127}},
{{4,32,216,360},{48,384,18,30},{288,36,300,180},{320,40,270,162}},
{{2,24,306,340},{180,15,330,297},{34,408,18,20},{396,33,150,135}},
{{9,34,192,396},{108,408,16,33},{306,81,330,160},{340,90,297,144}},
{{20,54,304,384},{160,432,38,48},{243,90,360,285},{405,150,216,171}},
{{17,39,312,432},{204,468,26,36},{351,153,360,260},{390,170,324,234}},
{{12,40,372,396},{144,480,31,33},{360,108,330,310},{400,120,297,279}},
{{6,48,324,540},{72,576,27,45},{432,54,450,270},{480,60,405,243}},
{{12,51,456,516},{144,612,38,43},{459,108,430,380},{510,120,387,342}},
{{8,53,348,600},{96,636,29,50},{477,72,500,290},{530,80,450,261}},
{{34,78,416,576},{272,624,52,72},{351,153,540,390},{585,255,324,234}},
{{24,80,496,528},{192,640,62,66},{360,108,495,465},{600,180,297,279}},
{{20,54,456,576},{240,648,38,48},{486,180,480,380},{540,200,432,342}},
{{17,55,288,648},{204,660,24,54},{495,153,540,240},{550,170,486,216}},
{{9,58,264,684},{108,696,22,57},{522,81,570,220},{580,90,513,198}},
{{50,414,486,538},{478,270,594,146},{542,178,218,550},{418,626,190,254}},
{{3,36,540,600},{264,22,590,531},{60,720,27,30},{708,59,220,198}},
{{8,64,432,720},{96,768,36,60},{576,72,600,360},{640,80,540,324}},
{{4,48,612,680},{360,30,660,594},{68,816,36,40},{792,66,300,270}},
{{24,102,608,688},{192,816,76,86},{459,108,645,570},{765,180,387,342}},
{{18,68,384,792},{216,816,32,66},{612,162,660,320},{680,180,594,288}},
{{205,515,551,625},{543,385,681,287},{627,301,345,623},{521,695,319,361}},
{{30,67,612,696},{360,804,51,58},{603,270,580,510},{670,300,522,459}},
{{16,106,464,800},{128,848,58,100},{477,72,750,435},{795,120,450,261}},
{{40,108,608,768},{320,864,76,96},{486,180,720,570},{810,300,432,342}},
{{34,110,384,864},{272,880,48,108},{495,153,810,360},{825,255,486,216}},
{{5,60,684,760},{76,912,45,50},{576,48,690,621},{828,69,480,432}},
{{17,76,456,876},{204,912,38,73},{684,153,730,380},{760,170,657,342}},
{{18,116,352,912},{144,928,44,114},{522,81,855,330},{870,135,513,198}},
{{20,54,646,816},{340,918,38,48},{405,150,792,627},{891,330,360,285}},
{{90,243,646,816},{432,160,792,627},{340,918,171,216},{891,330,384,304}},
{{6,48,540,900},{120,960,27,45},{352,44,885,531},{944,118,330,198}},
{{10,80,540,900},{120,960,45,75},{720,90,750,450},{800,100,675,405}},
{{34,78,624,864},{408,936,52,72},{702,306,720,520},{780,340,648,468}},
{{15,80,648,852},{180,960,54,71},{720,135,710,540},{800,150,639,486}},
{{24,80,744,792},{288,960,62,66},{720,216,660,620},{800,240,594,558}},
{{30,81,684,864},{360,972,57,72},{729,270,720,570},{810,300,648,513}}}
ArrayReshape[#^3,{4,4}]&/@Permutations[Range[16]]
but that's 20922789888000 items, so you stand no chance of doing it this way. Your best course of action would be to useFindInstance
as I described in a previous question of yours, but drop the diagonal constraints. $\endgroup$pwrs = ParallelTable[{n, PowersRepresentations[n, 4, 3]}, {n, 1, 4000}]; Select[pwrs, Length[Last[#]] >= 4 &]
-- edit: in fact there appear to be none all the way out to 18000 withLength[Last[#]] >= 8 &]
$\endgroup$