Pythagorean quadruples

Pythagorean quadruples is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

One form of   Pythagorean quadruples   is   (for positive integers   a,   b,   c,   and   d):

  a²   +   b²   +   c²     =     d²

An example:

  2²   +   3²   +   6²     =     7²

which is:

  4   +   9   +   36     =     49

Task

For positive integers up   2,200   (inclusive),   for all values of   a,   b,   and   c,
find   (and show here)   those values of   d   that   can't   be represented.

Show the values of   d   on one line of output   (optionally with a title).

Related task

•   Euler's sum of powers conjecture.

C

Five seconds on my Intel Linux box. <lang C>#include <stdio.h>

1. include <math.h>
2. include <string.h>

1. define N 2200

int main(int argc, char **argv){

  int a,b,c,d;
  int r[N+1];
  memset(r,0,sizeof(r));	// zero solution array
  for(a=1; a<=N; a++){
     for(b=a; b<=N; b++){

int aabb; if(a&1 && b&1) continue; // for positive odd a and b, no solution. aabb=a*a + b*b; for(c=b; c<=N; c++){ int aabbcc=aabb + c*c; d=(int)sqrt((float)aabbcc); if(aabbcc == d*d && d<=N) r[d]=1; // solution }

     }
  }
  for(a=1; a<=N; a++)
     if(!r[a]) printf("%d ",a);	// print non solution
  printf("\n");

}</lang>

Output:

$ clang -O3 foo.c -lm
$ ./a.out
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048 

Kotlin

This uses a similar approach to the REXX optimized version. It also takes advantage of a hint in the C entry that there is no solution if both a and b are odd (confirmed by Wikipedia article). Runs in about 7 seconds on my modest laptop which is more than 4 times faster than the brute force version would have been: <lang scala>// version 1.1.3

const val MAX = 2200 const val MAX2 = MAX * MAX - 1

fun main(args: Array<String>) {

   val found = BooleanArray(MAX + 1)       // all false by default
   val p2 = IntArray(MAX + 1) { it * it }  // pre-compute squares

   // compute all possible positive values of d * d - c * c and map them back to d
   val dc = mutableMapOf<Int, MutableList<Int>>()
   for (d in 1..MAX) {
       for (c in 1 until d) {
           val diff = p2[d] - p2[c]              
           val v = dc[diff]
           if (v == null)
               dc.put(diff, mutableListOf(d))
           else if (d !in v)
               v.add(d)
       }
   }
 
   for (a in 1..MAX) {
       for (b in 1..a) {
           if ((a and 1) != 0 && (b and 1) != 0) continue
           val sum = p2[a] + p2[b]
           if (sum > MAX2) continue
           val v = dc[sum]
           if (v != null) v.forEach { found[it] = true }
       }
   }
   println("The values of d <= $MAX which can't be represented:")
   for (i in 1..MAX) if (!found[i]) print("$i ")
   println()

}</lang>

Output:

The values of d <= 2200 which can't be represented:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048 

REXX

brute force

This version is a brute force algorithm, with some optimization (to save compute time)
which pre-computes some of the squares of the positive integers used in the search.
Integers too large for the precomputed values use an optimized integer square root. <lang rexx>/*REXX program computes/displays integers that aren't possible for: a² + b² + c² = d² */ parse arg hi . /*obtain optional argument from the CL.*/ if hi== | hi=="," then hi=2200; high=2*hi**2 /*Not specified? Then use the default.*/ @.=. /*array of integers to be squared. */ !.=. /* " " " squared. */

      do j=1  for high                          /*precompute possible squares (to max).*/
      jj=j*j;   !.jj=j;   if j<=hi  then @.j=jj /*define a square; D  value; squared # */
      end   /*j*/

d.=. /*array of possible solutions (D) */

      do       a=1  for hi-2                    /*go hunting for solutions to equation.*/
         do    b=a   to hi-1;   ab = @.a + @.b  /*calculate sum of  2  (A,B)   squares.*/
            do c=b   to hi;     abc= ab  + @.c  /*    "      "   "  3  (A,B,C)    "    */
            if abc<=high  then if !.abc==.  then iterate   /*Not a square? Then skip it*/
                          else do;   t=abc;   r=0;   q=1;       do while q<=t;   q=q*4
                                                                end   /*while q<=t*/
                                                           /*R:  will be the iSqrt(Z). */
                                 do while q>1;  q=q%4;  _=t-r-q;  r=r%2
                                 if _>=0  then do;  t=_;  r=r+q;  end
                                 end   /*while q>1*/
                               if r*r\=abc  then iterate
                               end
            s=!.abc;    d.s=                    /*define this D solution as being found*/
            end   /*c*/
         end      /*b*/
      end         /*a*/

say say 'Not possible positive integers for d ≤' hi " using equation: a² + b² + c² = d²" say $= /* [↓] find all the "not possibles". */

      do p=1  for hi;   if d.p==.  then $=$ p   /*Not possible? Then add it to the list*/
      end   /*p*/                               /* [↓]  display list of not-possibles. */

say substr($, 2) /*stick a fork in it, we're all done. */</lang>

output   when using the default input:

Not possible positive integers for   d ≤ 2200   using equation:  a² + b² + c²  =  d²

1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

optimized

This REXX version is an optimized version, it solves the formula:

  a²   +   b²     =     d²   -   c²

This REXX version is over thirty times faster then the previous version. <lang rexx>/*REXX program computes/displays integers that aren't possible for: a² + b² + c² = d² */ parse arg hi . /*obtain optional argument from the CL.*/ if hi== | hi=="," then hi=2200 /*Not specified? Then use the default.*/ high=hi * 3 /*D can be three times the HI (max).*/ @.=. /*array of integers (≤ hi) squared.*/ dbs.= /* " " differences between squares.*/

     do s=1  for high;  ss=s*s;  r.ss=s; @.s=ss /*precompute squares and square roots. */
     end  /*s*/

     do    c=1   for high;       cc=@.c         /*precompute possible differences.     */
        do d=c+1  to high;       dif=@.d - cc   /*process  D  squared; calc differences*/
        if wordpos(cc, dbs.dif)==0  then dbs.dif=dbs.dif cc     /*Not in list?  Add it.*/
        end   /*d*/
     end      /*c*/

d.=. /*array of possible solutions (D) */

     do     a=1  for hi-2                       /*go hunting for solutions to equation.*/
        do  b=a   to hi-1;         ab=@.a + @.b /*calculate sum of  2  (A,B)   squares.*/
        if dbs.ab==  then iterate             /*Not a difference?   Then ignore it.  */
           do n=1  for words(dbs.ab)            /*handle all ints that satisfy equation*/
           x=word(dbs.ab, n)                    /*get a C integer from the DBS.AB list.*/
           abc=ab + x                           /*add the  C²  integer  to  A²  +  B²  */
           _=r.abc                              /*retrieve the square root  of  C²     */
           d._=                                 /*mark the  D  integer as being found. */
           end   /*n*/
        end      /*b*/
     end         /*a*/

say say 'Not possible positive integers for d ≤' hi " using equation: a² + b² + c² = d²" say $= /* [↓] find all the "not possibles". */

      do p=1  for hi;   if d.p==.  then $=$ p   /*Not possible? Then add it to the list*/
      end   /*p*/                               /* [↓]  display list of not-possibles. */

say substr($, 2) /*stick a fork in it, we're all done. */</lang>

output   is the same as the 1^st REXX version.