How can I solve this system of 6 simple quadratic diophantine equations in four variables?

Question

I generated this text file containing tuples $(\sqrt{s}, \sqrt{t}, \sqrt{u}, s,t,u,t+u,t+u-s,t-s)$ where the six variables $s,t,u,t+u,t+u-s,t-s$ are all squares. An excerpt of the data set is:

153, 185, 672, 23409, 34225, 451584, 485809, 462400, 10816
264, 520, 975, 69696, 270400, 950625, 1221025, 1151329, 200704
264, 561, 952, 69696, 314721, 906304, 1221025, 1151329, 245025
306, 370, 1344, 93636, 136900, 1806336, 1943236, 1849600, 43264
357, 1325, 6960, 127449, 1755625, 48441600, 50197225, 50069776, 1628176
448, 952, 495, 200704, 906304, 245025, 1151329, 950625, 705600
459, 555, 2016, 210681, 308025, 4064256, 4372281, 4161600, 97344

My task is to iterate each line/row of this file and to try finding a solution $0<w<x<y<z$ of the following system of six diophantine equation:

• $x^2-w^2=s\qquad z^2-w^2=t+u$
• $y^2-w^2=t\qquad z^2-x^2=t+u-s$
• $z^2-y^2=u\qquad y^2-x^2=t-s$

What I tried so far is running the following script with more relaxed conditions:

ClearAll["Global`*"]
arr = Import[
   "C:/Users/esultano/git/pythagorean/pythagorean_stu_Arty_.txt", 
   "CSV", "HeaderLines" -> 0];
f[i_] := Part[arr[[i]], 1 ;; 3];
len = Length[arr];
For[i = 1, i < len, i++, {
  triplet = f[i];
  s = triplet[[1]];
  t = triplet[[2]];
  u = triplet[[3]];
  (*Print[FactorInteger[s]];*)
  ins = FindInstance[
    x*x - w*w == s*s && y*y - w*w == t*t && w != 0, {w, x, y}, 
    Integers];
  If[Length[ins] > 
    0, {ins = 
     Join[First[ins], {"z" -> N[Sqrt[u^2 + First[ins][[3, 2]]^2], 20],
        "s" -> s, "t" -> t, "u" -> u}]; Print[ins]}, Continue];
  ins = FindInstance[
    y*y - w*w == t*t && z*z - y*y == u*u && w != 0, {w, y, z}, 
    Integers];
  If[Length[ins] > 
    0, {ins = 
     Join[First[ins], {"x" -> N[Sqrt[s^2 + First[ins][[1, 2]]^2], 20],
        "s" -> s, "t" -> t, "u" -> u}]; Print[ins]}, Continue];
  }
 ]

It gives me "almost" solutions, for example $(w,x,y,z)=(40579,58565,65221,196605.2940)$. It seems not to be possible to combine the conditions such for example:

ins = FindInstance[x*x-w*w==s*s && y*y-w*w==t*t && z*z-y*y==u*u && w > 0, {w, x, y, z}, Integers]

Is there a way to let Mathematica solve such a system of diophantine equations?

The background is a quite inteersting problem, namely a search for four squares as described here at this MSE post - an extended version of Mengoli's Six Square Problem. It would be great if Mathematica can solve the above shown system of these six diophantine equations.

Answer 1

ClearAll["Global`*"]

arr = Cases[
   Import[
    "/Users/roberthanlon/Downloads/data.txt", "CSV"],
   {_Integer ..}, 1];

Length@arr

(* 27060 *)

Verifying that every line of arr matches the template {Sqrt[s], Sqrt[t], Sqrt[u], s, t, u, t+u, t+u-s, t-s}

And @@ ((repl = Thread[{s, t, u} -> #[[4 ;; 6]]];
     # == ({Sqrt[s], Sqrt[t], Sqrt[u], s, t, u, t + u, t + u - s, t - s} /.
              repl)) & /@ arr)

(* True *)

Since there are more equations than variables, the Solve option MaxExtraConditions must be used.

sol[s_Integer?Positive, t_Integer?Positive, u_Integer?Positive] :=
 Solve[{
   x^2 - w^2 == s, z^2 - w^2 == t + u, y^2 - w^2 == t,
   z^2 - x^2 == t + u - s, z^2 - y^2 == u, y^2 - x^2 == t - s,
   0 < w < x < y < z}, {w, x, y, z}, Integers,
  MaxExtraConditions -> All]

There are no solutions for the first 20 lines of arr

at = AbsoluteTiming[
  arr2 = Select[arr[[1 ;; 20]], Length[sol @@ #[[4 ;; 6]]] > 0 &]]

(* {2.05478, {}} *)

It would take about 50 minutes to run through the data (assuming that nothing else is competing for processing time).

Length[arr]*(at/20)/60

(* 46.3353 *)