4X4 matrix of distinct cubed integers with identical row and column sums
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}}}
All credit to @Roman who deserves the accept for discovering 7095816 was the target number. With this number pinned down I tried to solve this myself with Mathematica, but it was still too slow:
variables = Array[x, 16];
matrix = ArrayReshape[variables, {4, 4}];
target = 7095816;
uniqueConstraint = (And @@ (Unequal @@@ Subsets[variables, {2}]));
rangeConstraint = (And @@ (1 <= # < 200 & /@ variables));
magicConstraint = And[
And @@ (# == target & /@ Total[Transpose[matrix^3]]),
And @@ (# == target & /@ Total[matrix^3])];
FindInstance[magicConstraint && uniqueConstraint &&
rangeConstraint, variables, PositiveIntegers]
However, I was able to solve it in around two minutes using the Z3 solver for Python. I figured since you have posted a handful of questions previously about these sorts of magic square problems, that you might find this useful should you wish to experiment further:
from z3 import *
s = Solver()
# the target sum
target = 7095816
# NxN matrix
n = 4
vs = [ Int("x_%s" % i) for i in range(n*n) ]
# variables are 1 to 200
cnumber = [ And(1 <= v, v <= 200) for v in vs ]
s.add(cnumber)
# all distinct
cdistinct = Distinct(*vs)
s.add(cdistinct)
def cubesum(vlist):
return sum(v*v*v for v in vlist)
# row totals equal column totals
for i in range(n):
s.add(cubesum(vs[i*n:(i+1)*n]) == target)
s.add(cubesum(vs[i::n]) == target)
if s.check() == sat:
m = s.model()
r = [m.evaluate(v) for v in vs]
print("sum: {}, vars: {}".format(target,r))
else:
print("failed")
This produced the result:
sum: 7095816, vars: [144, 18, 150, 90, 160, 20, 135, 81, 2, 16, 108, 180, 24, 192, 9, 15]
... and if you look closely this solution is a permuation of @Roman's matrix.