4X4 matrix of distinct cubed integers with identical row and column sums

Question

I wrote the following code:

ClearAll[findMat]
findMat = 
  Module[{csums = ConstantArray[#, 4], 
     mats = Select[DuplicateFreeQ[Join @@ #] &]@
       Subsets[Select[DuplicateFreeQ]@
         IntegerPartitions[#, {4}, Range[Ceiling[#^(1/3)]]^3], {4}]}, 
    Join @@ ((Select[Total[#] == csums &]@
          Tuples[{{#[[1]]}, Permutations[#[[2]]], 
            Permutations[#[[3]]], 
            Permutations[#[[4]]]}]) & /@ mats)] &;
solsB = ParallelTable[
   findMat2[n] /. {} -> Nothing, {n, 0, 10^5}];
{Total[#[[1, 1]]], MatrixForm /@ #} & /@ solsB

This code looks for a 4X4 matrix of distinct cubed integers with identical row and column sums. But Mathematica gave me an error about insufficient memory. How can I solve this problem? Or did I make mistake in my code?

Your code is very terse and compact. It may help explain what you're doing if you add some comments. What range of numbers can go in the matrix? If it's 1 to 16, then you need to test all of 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 use FindInstance as I described in a previous question of yours, but drop the diagonal constraints. — flinty, Commented Jul 24, 2020 at 11:23
I commented the code there - you just need to delete two lines - the only other thing is you need to know the total in advance, but that could be removed if I introduce another constraint that all totals are equal mathematica.stackexchange.com/a/226351/72682 — flinty, Commented Jul 24, 2020 at 11:28
@flinty So I need to delete the stuf under: * Both diagonals of the matrix *? — Jan Eerland, Commented Jul 24, 2020 at 11:31
The trouble is that there are very few small integers that can be written as a sum of four cubes in 8 ways which is a necessary but not sufficient requirement: 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 with Length[Last[#]] >= 8 &] — flinty, Commented Jul 24, 2020 at 11:57

Roman · Accepted Answer · 2020-08-10 17:17:47Z

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}}}

+1, you put a lot of effort into this. How long was it running for? You should also add -std=c99 and -lm to the compile options, because otherwise it will complain about libmath linkage. — flinty, Commented Jul 31, 2020 at 13:33
@flinty thanks for the options, they weren't necessary on macOS so I didn't bother, but Linux may be different. The code took several hours to run on my 2.4 GHz Intel Core i9 (MacBook Pro). — Roman, Commented Jul 31, 2020 at 13:51

flinty · Accepted Answer · 2020-07-31 18:14:57Z

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.

Stack Exchange Network

4X4 matrix of distinct cubed integers with identical row and column sums

2 Answers 2

C code

Update

Your Answer

Not the answer you're looking for? Browse other questions tagged
calculus-and-analysis
matrix
parallelization
code-review
permutation
or ask your own question.

Linked

Hot Network Questions

4X4 matrix of distinct cubed integers with identical row and column sums

2 Answers 2

C code

Update

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged calculus-and-analysismatrixparallelizationcode-reviewpermutation or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
calculus-and-analysis
matrix
parallelization
code-review
permutation
or ask your own question.