Pythagorean quadruples

From Rosetta Code
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):


  a2   +   b2   +   c2     =     d2


An example:

  22   +   32   +   62     =     72
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 tasks



ALGOL 68

As with the optimised REXX solution, we find the values of d for which there are no a^2 + b^2 = d^2 - c^2.

BEGIN
# find values of d where d^2 =/= a^2 + b^2 + c^2 for any integers a, b, c #
# where d in [1..2200], a, b, c =/= 0 #
# max number to check #
INT max number = 2200;
INT max square = max number * max number;
# table of numbers that can be the sum of two squares #
[ 1 : max square ]BOOL sum of two squares; FOR n TO max square DO sum of two squares[ n ] := FALSE OD;
FOR a TO max number DO
INT a2 = a * a;
FOR b FROM a TO max number WHILE INT sum2 = ( b * b ) + a2;
sum2 <= max square DO
sum of two squares[ sum2 ] := TRUE
OD
OD;
# now find d such that d^2 - c^2 is in sum of two squares #
[ 1 : max number ]BOOL solution; FOR n TO max number DO solution[ n ] := FALSE OD;
FOR d TO max number DO
INT d2 = d * d;
FOR c TO d - 1 WHILE NOT solution[ d ] DO
INT diff2 = d2 - ( c * c );
IF sum of two squares[ diff2 ] THEN
solution[ d ] := TRUE
FI
OD
OD;
# print the numbers whose squares are not the sum of three squares #
FOR d TO max number DO
IF NOT solution[ d ] THEN
print( ( " ", whole( d, 0 ) ) )
FI
OD;
print( ( newline ) )
END
Output:
 1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

C

Version 1

Five seconds on my Intel Linux box.

#include <stdio.h>
#include <math.h>
#include <string.h>
 
#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");
}
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 

Version 2 (much faster)

Translation of second version of FreeBASIC entry. Runs in about 0.15 seconds.

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
 
#define N 2200
#define N2 2200 * 2200 * 2
 
int main(int argc, char **argv) {
int a, b, c, d, a2, s = 3, s1, s2;
int r[N + 1];
memset(r, 0, sizeof(r));
int *ab = calloc(N2 + 1, sizeof(int)); // allocate on heap, zero filled
 
for (a = 1; a <= N; a++) {
a2 = a * a;
for (b = a; b <= N; b++) ab[a2 + b * b] = 1;
}
 
for (c = 1; c <= N; c++) {
s1 = s;
s += 2;
s2 = s;
for (d = c + 1; d <= N; d++) {
if (ab[s1]) r[d] = 1;
s1 += s2;
s2 += 2;
}
}
 
for (d = 1; d <= N; d++) {
if (!r[d]) printf("%d ", d);
}
printf("\n");
free(ab);
return 0;
}
Output:
Same as first version.

FreeBASIC

From the Wikipedia page

Alternate parametrization, second version both A and B even.Time just less then 0.9 second on a AMD Athlon II X4 645 3.34GHz win7 64bit. Program uses one core. When the limit is set to 576 (abs. minimum for 2200), the time is about 0.85 sec.

' version 16-07-2017
' compile with: fbc -s console
 
#Define max 2200
 
Dim As ULong l, m, n, l2, l2m2,list_1(), list_2()
Dim As ULong limit = max * 4 \ 15
Dim As ULong max2 = limit * limit * 2
ReDim As ULong list_1(max2), list_2(max2)
 
' prime sieve, list_2(l) contains a 0 if l = prime
For l = 4 To max2 Step 2
list_1(l) = 1
Next
For l = 3 To max2 Step 2
If list_1(l) = 0 Then
For m = l * l To max2 Step l * 2
list_1(m) = 1
Next
End If
Next
 
' we do not need a and b (a and b are even, l = a \ 2, m = b \ 2)
' we only need to find d
For l = 1 To limit
l2 = l * l
For m = l To limit
l2m2 = l2 + m * m
list_2(l2m2 +1) = 1
' if l2m2 is a prime, no other factors exits,
If list_1(l2m2) = 0 Then Continue For
' find possible factors of l2m2
For n = 2 To Fix(Sqr(l2m2 -1))
If l2m2 Mod n = 0 Then
' set list_2(x) to 1 if solution is found
list_2(l2m2 \ n + n) = 1
End If
Next
Next
Next
 
For l = 1 To max
If list_2(l) = 0 Then Print l; " ";
Next
Print
 
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

Brute force

Based on the second REXX version: A^2 + B^2 = D^2 - C^2. Faster then the first version, about 0.2 second

' version 21-07-2017
' compile with: fbc -s console
 
#Define n 2200
 
Dim As ULong s = 3, s1, s2, x, x2, y
ReDim As ULong l(n), l_add(n * n * 2)
 
For x = 1 To n
x2 = x * x
For y = x To n
l_add(x2 + y * y) = 1
Next
Next
 
For x = 1 To n
s1 = s
s += 2
s2 = s
For y = x +1 To n
If l_add(s1) = 1 Then l(y) = 1
s1 += s2
s2 += 2
Next
Next
 
For x = 1 To n
If l(x) = 0 Then Print Str(x); " ";
Next
Print
 
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

Kotlin

Version 1

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:

// 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()
}
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 

Version 2 (much faster)

This is a translation of the second FreeBASIC version and runs in about the same time (0.2 seconds).

One thing I've noticed about the resulting sequence is that it appears to be an interleaving of the two series 2 ^ n and 5 * (2 ^ n) for n >= 0 though whether it's possible to prove this mathematically I don't know.

// version 1.1.3
 
const val MAX = 2200
const val MAX2 = MAX * MAX * 2
 
fun main(args: Array<String>) {
val found = BooleanArray(MAX + 1) // all false by default
val a2b2 = BooleanArray(MAX2 + 1) // ditto
var s = 3
 
for (a in 1..MAX) {
val a2 = a * a
for (b in a..MAX) a2b2[a2 + b * b] = true
}
 
for (c in 1..MAX) {
var s1 = s
s += 2
var s2 = s
for (d in (c + 1)..MAX) {
if (a2b2[s1]) found[d] = true
s1 += s2
s2 += 2
}
}
 
println("The values of d <= $MAX which can't be represented:")
for (d in 1..MAX) if (!found[d]) print("$d ")
println()
}
Output:
Same as Version 1.

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.

/*REXX pgm computes/shows (integers),  D  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(t). */
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. */
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:

  a2   +   b2     =     d2   -   c2

This REXX version is over thirty times faster then the previous version.

/*REXX pgm computes/shows (integers),  D  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. */
output   is the same as the 1st REXX version.


zkl

Translation of: ALGOL 68
# find values of d where d^2 =/= a^2 + b^2 + c^2 for any integers a, b, c #
# where d in [1..2200], a, b, c =/= 0 #
# max number to check #
const max_number = 2200;
const max_square = max_number * max_number;
# table of numbers that can be the sum of two squares #
sum_of_two_squares:=Data(max_square+1,Int).fill(0); # 4 meg byte array
foreach a in ([1..max_number]){
a2 := a * a;
foreach b in ([a..max_number]){
sum2 := ( b * b ) + a2;
if(sum2 <= max_square) sum_of_two_squares[ sum2 ] = True; # True-->1
}
}
# now find d such that d^2 - c^2 is in sum of two squares #
solution:=Data(max_number+1,Int).fill(0); # another byte array
foreach d in ([1..max_number]){
d2 := d * d;
foreach c in ([1..d-1]){
diff2 := d2 - ( c * c );
if(sum_of_two_squares[ diff2 ]){ solution[ d ] = True; break; }
}
}
# print the numbers whose squares are not the sum of three squares #
foreach d in ([1..max_number]){
if(not solution[ d ]) print(d, " ");
}
println();
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048