Pythagorean quadruples

From Rosetta Code
Task
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[edit]

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[edit]

Version 1[edit]

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)[edit]

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[edit]

From the Wikipedia page[edit]

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.

' 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
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

Brute force[edit]

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

' 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
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

Go[edit]

Translation of: FreeBASIC
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()
}
Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048 

jq[edit]

The following is a direct solution but with some obvious optimizations. Its main value may be to illustrate how looping with breaks can be accomplished in jq without `foreach`. Notice also how `first/1` is used in `is_pythagorean_quad/0` to avoid unnecessary computation.

# Emit a proof that the input is a pythagorean quad, or else false
def is_pythagorean_quad:
. as $d
| (.*.) as $d2
| first(
label $continue_a | range(1; $d) | . as $a | (.*.) as $a2
| if 3*$a2 > $d2 then break $continue_a else . end
| label $continue_b | range($a; $d) | . as $b | (.*.) as $b2
| if $a2 + 2 * $b2 > $d2 then break $continue_b else . end
| (($d2-($a2+$b2)) | sqrt) as $c
| if ($c | floor) == $c then [$a, $b, $c] else empty end )
// false;
 
# The specific task:
 
[range(1; 2201) | select( is_pythagorean_quad | not )] | join(" ")

Invocation and Output

jq -r -n -f program.jq
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048

Kotlin[edit]

Version 1[edit]

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)[edit]

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.

Perl 6[edit]

Works with: Rakudo version 2017.08
my \N = 2200;
my @sq = (0 .. N)»²;
my @not = False xx N;
@not[0] = True;
 
for 1 .. N -> $d {
my $last = 0;
for $d ... ($d/3).ceiling -> $a {
for 1 .. ($a/2).ceiling -> $b {
last if (my $ab = @sq[$a] + @sq[$b]) > @sq[$d];
if (@sq[$d] - $ab).sqrt.narrow ~~ Int {
@not[$d] = True;
$last = 1;
last
}
}
last if $last;
}
}
 
say @not.grep( *.not, :k );
Output:
(1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048)

REXX[edit]

brute force[edit]

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[edit]

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; [email protected].c /*precompute possible differences. */
do d=c+1 to high; [email protected].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; [email protected].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[edit]

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