Pythagorean quadruples
You are encouraged to solve this task according to the task description, using any language you may know.
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, c, and d,
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. <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
Version 1
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
Version 2 (much faster)
Translation of second version of FreeBASIC entry. Runs in about 0.15 seconds. <lang C>#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;
}</lang>
- Output:
Same as first version.
FreeBASIC
From the Wikipedia page
Alternate parametrization, second version both A and B even.Time just less then 0.7 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. <lang freebasic>' version 12-08-2017 ' compile with: fbc -s console
- Define max 2200
Dim As UInteger l, m, n, l2, l2m2 Dim As UInteger limit = max * 4 \ 15 Dim As UInteger max2 = limit * limit * 2 ReDim As Ubyte list_1(max2), list_2(max2 +1)
' 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
' if l2m2 is odd, we need only to check the odd divisors
For n = 2 + (l2m2 And 1) To Fix(Sqr(l2m2 -1)) Step 1 + (l2m2 And 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</lang>
- 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 <lang freebasic>' version 14-08-2017 ' compile with: fbc -s console
- Define n 2200
Dim As UInteger s = 3, s1, s2, x, x2, y ReDim As Ubyte 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</lang>
- Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048
Go
<lang go>package main
import "fmt"
const (
N = 2200 N2 = N * N * 2
)
func main() {
s := 3 var s1, s2 int var r [N + 1]bool var ab [N2 + 1]bool
for a := 1; a <= N; a++ {
a2 := a * a
for b := a; b <= N; b++ {
ab[a2 + b * b] = true
}
}
for c := 1; c <= N; c++ {
s1 = s
s += 2
s2 = s
for d := c + 1; d <= N; d++ {
if ab[s1] {
r[d] = true
}
s1 += s2
s2 += 2
}
}
for d := 1; d <= N; d++ {
if !r[d] {
fmt.Printf("%d ", d)
}
}
fmt.Println()
}</lang>
- 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: <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
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. <lang scala>// 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()
}</lang>
- 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.
<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(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. */</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