Simultaneously Square Sums of Mostly Squares

Background

This problem is one of a family of related problems:

a²+b²+c²=d²

The well-known Pythagorean Quadruples problem.

a²+b²+c=d² for all (a,b,c) permutations

This is what you are here for!

a²+b+c=d² for all (a,b,c) permutations

A surprisingly simple Asia-Pacific Mathematical Olympiad 2011 question as solved in Perfect Squares are Impossible - APMO 2011.

a+b+c=d²

Wat? That's not a problem.

Exposition

So we're looking for triples of positive integers a,b,c where all of three permutations of the not-quite-sums-of-squares a+b²+c², a²+b+c², and a²+b²+c are perfect squares.

The fully generic problem is the most interesting one, but there will be interesting - or boring, depending on your point of view - sub-families of results where two or even three of the variables are equal.

If there was nothing special about the distribution of solutions in ℤ³, then the "volume" of the search space for each of these classes is incomparable. The a=b=c case is just a line, the a=b and b=c cases are planes, and the generic case is the full rest of the solid region. One might therefore expect almost none of the solutions to be the special cases, but those special properties might enable them to have an unexpected density to their solutions.

Of course, the solutions where these planes intersect is precisely the set of solutions along the diagonal; and of course, the solutions on the planes are precisely the intersection of the planes with the solutions to the generic version of the problem. There is only one problem being solved, we're just slicing it different ways.

a=b=c

Note that the a=b=c case is a simple quadratic diophantine equation 2·a²+a=d², expressable as a pell equation (4·a+1)²-2·(2·d)²=1. The infinite family of solutions for this equation is known, they are derived from the subset of the solutions to the classic Pell's Equation x²-2·y²=1, with x restricted to numbers of the form 4a+1, and y even. This restriction gives them more in common with the approximants to √32, than those of just √2. This gives them a known density, they grow, and thus thin out, exponentially. They are marked up thusly in the results table.

a=b xor b=c

The a=b and b=c cases (by symmetry, these are the same case, but I've ordeded a,b, and c, so they superficially fall into two subfamilies: one-small-two-big, and two-small-one-big) also have a well-studied algebraic property, based on a similar substitution as performed above. a²+a+c²=d² is precisely 4·a²+4·a+1+4·c^2=4·d^2+1, so (2·a+1)^2+(2·c)^2=(2·d)^2+1, so this is an instance of the "x²+y²-z²=1" problem, the hyperboloid of one sheet, with solutions as documented out there in internet land.

Of course, our solutions also need to satisfy 2·a²+c=f², so are a very small subset of the solutions to that simpler problem. Some proper mathematician needs to have a look at this, I don't have the chops.

a≠b≠c

No clue - if you can extract anything from the known solutions, or hit this with actual maths, drop me a note. The solutions a≠b≠c are higlighted thus in the table of results.

One thing's certain - this is tractible, and the right mind will probably solve it easily. In fact, it seems Diophantus would have solved it, as he solved a harder problem, Diophantus Problem III.19 (page 8).

Results

Code

Multi-threaded C

All-round guru Robert Gerbiecz from the Mersenne Forum completely reorganised Inz' C code to utilise a small-prime sieve to filter out large swathes of impossible search space. It seems to provide a ~5x improvement over the original C (still being tested).

Inz converted the GP code into a multithreaded C program at something like 100x speedup, or something absurd like that. Compile it with:

gcc -Wall -Wextra -O3 -o sssms2 sssms2.c -lm -lpthread

Pari/GP Script

Here's a script that will search the range 1≤c≤b≤a≤5000. Those bounds being the first two numbers in the outermost 'for' command:

f2=[\
        19,1,0,0,17,0,0,0;\
        1,0,0,0,0,0,0,0;\
        0,0,0,136,0,0,0,136;\
        0,0,136,68,0,0,136,68;\
        17,0,0,0,17,0,0,0;\
        0,0,0,0,0,0,0,0;\
        0,0,0,136,0,0,0,136;\
        0,0,136,68,0,0,136,196\
  ]
  for(a=1,5000,\
      a8=a%8;\
      for(b=1,a,\
          b8=b%8;\
          f3=f2[1+a8,1+b8];\
          if(f3,\
             a2b2=a*a+b*b;\
             mrc=floor(sqrt(a2b2+b));\
             for(rc=sqrtint(a2b2-1)+1,mrc,\
                 c=rc*rc-a2b2;\
                 if(bittest(f3,c%8)&&\
                    issquare(a*a+b+c*c)&&\
                    issquare(a+b*b+c*c),\
                    print(c" "b" "a)\
                 )\
              )\
          )\
    )\
)

As one-line commands:

f2=[19,1,0,0,17,0,0,0;1,0,0,0,0,0,0,0;0,0,0,136,0,0,0,136;0,0,136,68,0,0,136,68;17,0,0,0,17,0,0,0;0,0,0,0,0,0,0,0;0,0,0,136,0,0,0,136;0,0,136,68,0,0,136,196]
dorange(lo,hi)=for(a=lo,hi,a8=a%8;for(b=1,a,b8=b%8;f3=f2[1+a8,1+b8];if(f3,a2b2=a*a+b*b;mrc=floor(sqrt(a2b2+b));for(rc=sqrtint(a2b2-1)+1,mrc,c=rc*rc-a2b2;if(bittest(f3,c%8)&&issquare(a*a+b+c*c)&&issquare(a+b*b+c*c),print(c" "b" "a))))))

which, if saved as sssms.gp, can be invoked thus:

echo 'dorange(1,5000)' | gp sssms.gp

or even:

echo -e '#\ndorange(1,5000)' | gp sssms.gp

I ran with the '#' directive in GP before those lines to turn on timing, so I could have a comparison of how quickly different versions of the code ran on different boxes. With my little box (a 1.5GHz Celeron, yup, in 2021) the above 1-5000 took time = 13,949 ms. Bands of 10000 at higher ranges are starting to take many hours. There's an argument for either using C or some more mathematical smarts to solve this.

Improvements Over Brute Force

Restrictions Modulo 8

The 8×8 array is a set of bitmasks of valid residues mod 8 that the third variable can have given the residues mod 8 that the other 2 variables have. This is particularly useful as squares can only be 0,1,and 4 (mod 8), and to get the 3 sums to all be simultaniously from those valid values discards many b and c values. Only 36 of the 8×8×8=512 possibilities are permitted. It even throws away 1/8 of the a values - those 5 mod 8, hence the immediate rejection of such values in the C code. There's no need to immediately reject b==5 (mod 8) separately, that's part of the test on f2[a&7][b&7].

Restrictions Modulo 3

None implemented so far, but it appears not all (a,b,c) triplets mod 3 are admissible either (not surprising, -2 isn't a square mod 3 so you'd expect only (2/3)³=8/27 of the candidates to survive such a test:

$ perl -pe 's/(\d+)\s+(\d+)\s+(\d+).*/($1%3)."\t".($2%3)."\t".($3%3)/e' < squaresums.list  | sort | uniq -c
     18 0	0	0
     42 0	0	1
      9 0	1	0
     20 1	0	0
     17 1	1	1
     19 1	1	2
      2 1	2	1
      8 1	2	2
      6 2	1	1
      5 2	1	2
     18 2	2	1
     43 2	2	2

So (000), (100)*, (111), (112)*, (122)*, and (222) are fine.

Therefore a==0 (mod 3) lets you reject b==2 (mod 3), and by symmetry, a==2 (mod 3) lets you reject b==0 (mod 3). a==1 (mod 3) brings no immediate gains.

Fin

Thanks to inz, boru, janrinok, and chromas from SoylentNews for letting me run the scripts/program on way bigger boxes than anything I have access to at home. And extra thanks to inz for converting the GP into C.