Pythagorean quadruples
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 task
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. <lang algol68>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</lang>
- Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048
C
Five seconds on my Intel Linux box. <lang C>#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");
}</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 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(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:
- a2 + b2 = d2 - c2
This REXX version is over thirty times faster then the previous version. <lang rexx>/*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. */</lang>
- output is the same as the 1st REXX version.
zkl
<lang zkl># 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();</lang>
- Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048