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.
D
<lang D>import std.bitmanip : BitArray; import std.stdio;
enum N = 2_200; enum N2 = 2*N*N;
void main() {
BitArray found; found.length = N+1;
BitArray aabb; aabb.length = N2+1;
uint s=3;
for (uint a=1; a<=N; ++a) {
uint aa = a*a;
for (uint b=1; b<N; ++b) {
aabb[aa + b*b] = true;
}
}
for (uint c=1; c<=N; ++c) {
uint s1 = s;
s += 2;
uint s2 = s;
for (uint d=c+1; d<=N; ++d) {
if (aabb[s1]) {
found[d] = true;
}
s1 += s2;
s2 += 2;
}
}
writeln("The values of d <= ", N, " which can't be represented:");
for (uint d=1; d<=N; ++d) {
if (!found[d]) {
write(d, ' ');
}
}
writeln;
}</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
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
Haskell
<lang Haskell>module PythagoreanQuadruple where
powersOfTwo = 1 : map (*2) powersOfTwo
unrepresentable = merge powersOfTwo $ map (*5) powersOfTwo
merge (x:xs) (y:ys) | x < y = x : merge xs (y:ys) merge (x:xs) (y:ys) = y : merge (x:xs) ys
main = do
putStrLn "The values of d <= 2200 which can't be represented." print $ takeWhile (<= 2200) unrepresentable
</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]
Java
<lang java>public class PythagoreanQuadruple {
static final int MAX = 2200; static final int MAX2 = MAX * MAX * 2;
public static void main(String[] args) {
boolean[] found = new boolean[MAX + 1]; // all false by default
boolean[] a2b2 = new boolean[MAX2 + 1]; // ditto
int s = 3;
for (int a = 1; a <= MAX; a++) {
int a2 = a * a;
for (int b = a; b <= MAX; b++) a2b2[a2 + b * b] = true;
}
for (int c = 1; c <= MAX; c++) {
int s1 = s;
s += 2;
int s2 = s;
for (int d = c + 1; d <= MAX; d++) {
if (a2b2[s1]) found[d] = true;
s1 += s2;
s2 += 2;
}
}
System.out.printf("The values of d <= %d which can't be represented:\n", MAX);
for (int d = 1; d <= MAX; d++) {
if (!found[d]) System.out.printf("%d ", d);
}
System.out.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
jq
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. <lang jq># 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(" ")</lang>
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
Julia
<lang julia>function quadruples(N::Int=2200)
r = falses(N) ab = falses(2N ^ 2)
for a in 1:N, b in a:N
ab[a ^ 2 + b ^ 2] = true
end
s = 3
for c in 1:N
s1, s, s2 = s, s + 2, s + 2
for d in c+1:N
if ab[s1] r[d] = true end
s1 += s2
s2 += 2
end
end
return find(.! r)
end
println("Pythagorean quadruples up to 2200: ", join(quadruples(), ", "))</lang>
- Output:
Pythagorean quadruples up to 2200: 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.
Perl 6
<lang perl6>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 );</lang>
- Output:
(1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048)
PicoLisp
<lang PicoLisp>(de quadruples (N)
(let (AB NIL S 3 R)
(for A N
(for (B A (>= N B) (inc B))
(idx
'AB
(+ (* A A) (* B B))
T ) ) )
(for C N
(let (S1 S S2)
(inc 'S 2)
(setq S2 S)
(for (D (+ C 1) (>= N D) (inc D))
(and (idx 'AB S1) (idx 'R D T))
(inc 'S1 S2)
(inc 'S2 2) ) ) )
(make
(for A N
(or (idx 'R A) (link A)) ) ) ) )
(println (quadruples 2200))</lang>
- Output:
(1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048)
Python
<lang Python>def quad(top=2200):
r = [False] * top
ab = [False] * (top * 2)**2
for a in range(1, top):
for b in range(a, top):
ab[a * a + b * b] = True
s = 3
for c in range(1, top):
s1, s, s2 = s, s + 2, s + 2
for d in range(c + 1, top):
if ab[s1]:
r[d] = True
s1 += s2
s2 += 2
return [i for i, val in enumerate(r) if not val and i]
if __name__ == '__main__':
n = 2200
print(f"Those values of d in 1..{n} that can't be represented: {quad(n)}")</lang>
- Output:
Those values of d in 1..2200 that 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]
Racket
<lang racket>#lang racket
(require data/bit-vector)
(define (quadruples top)
(define top+1 (add1 top)) (define 1..top (in-range 1 top+1)) (define r (make-bit-vector top+1)) (define ab (make-bit-vector (add1 (sqr (* top 2))))) (for* ((a 1..top) (b (in-range a top+1))) (bit-vector-set! ab (+ (sqr a) (sqr b)) #t))
(for/fold ((s 3))
((c 1..top))
(for/fold ((s1 s) (s2 (+ s 2)))
((d (in-range (add1 c) top+1)))
(when (bit-vector-ref ab s1)
(bit-vector-set! r d #t))
(values (+ s1 s2) (+ s2 2)))
(+ 2 s))
(for/list ((i (in-naturals 1)) (v (in-bit-vector r 1)) #:unless v) i))
(define (report n)
(printf "Those values of d in 1..~a that can't be represented: ~a~%" n (quadruples n)))
(report 2200)</lang>
- Output:
Those values of d in 1..2200 that 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(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 /*(DO iterator C)*/
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.
Ring
<lang ring># Project : Pythagorean quadruples
- Date : 2017/12/17
- Author : Gal Zsolt (~ CalmoSoft ~)
- Email : <calmosoft@gmail.com>
limit = 2200 pq = list(limit) for n = 1 to limit
for m = 1 to limit
for p = 1 to limit
for x = 1 to limit
if pow(x,2) = pow(n,2) + pow(m,2) + pow(p,2)
pq[x] = 1
ok
next
next
next
next pqstr = "" for d = 1 to limit
if pq[d] = 0
pqstr = pqstr + d + " "
ok
next see pqstr + nl </lang>
- Output:
1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048
Scala
- Output:
Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).
<lang Scala>object PythagoreanQuadruple extends App {
val MAX = 2200
val MAX2: Int = MAX * MAX * 2
val found = Array.ofDim[Boolean](MAX + 1)
val a2b2 = Array.ofDim[Boolean](MAX2 + 1)
var s = 3
for (a <- 1 to MAX) {
val a2 = a * a
for (b <- a to MAX) a2b2(a2 + b * b) = true }
for (c <- 1 to MAX) {
var s1 = s
s += 2
var s2 = s
for (d <- (c + 1) to MAX) {
if (a2b2(s1)) found(d) = true
s1 += s2
s2 += 2
}
}
println(f"The values of d <= ${MAX}%d which can't be represented:")
val notRepresented = (1 to MAX).filterNot(d => found(d) )
println(notRepresented.mkString(" "))
}</lang>
Sidef
<lang ruby>require('ntheory')
- Finds all solutions (a,b) such that: a^2 + b^2 = n^2
func sum_of_two_squares(n) is cached {
n == 0 && return 0, 0
var prod1 = 1 var prod2 = 1
var prime_powers = []
for p,e in (n.factor_exp) {
if (p % 4 == 3) { # p = 3 (mod 4)
e.is_even || return [] # power must be even
prod2 *= p**(e >> 1)
}
elsif (p == 2) { # p = 2
if (e.is_even) { # power is even
prod2 *= p**(e >> 1)
}
else { # power is odd
prod1 *= p
prod2 *= p**((e - 1) >> 1)
prime_powers.append([p, 1])
}
}
else { # p = 1 (mod 4)
prod1 *= p**e
prime_powers.append([p, e])
}
}
prod1 == 1 && return prod2, 0 prod1 == 2 && return prod2, prod2
# All the solutions to the congruence: x^2 = -1 (mod prod1)
var square_roots = gather {
gather {
for p,e in (prime_powers) {
var pp = p**e
var r = %S<ntheory>.sqrtmod("#{pp - 1}", "#{pp}")
take([[r, pp], [pp - r, pp]])
}
}.cartesian { |*a|
take(Math.chinese(a...))
}
}
var solutions = []
for r in (square_roots) {
var s = r
var q = prod1
while (s*s > prod1) {
(s, q) = (q % s, s)
}
solutions.append([prod2 * s, prod2 * (q % s)]) }
for p,e in (prime_powers) {
for (var i = e%2; i < e; i += 2) {
var sq = p**((e - i) >> 1)
var pp = p**(e - i)
solutions += (
__FUNC__(prod1 / pp).map { |pair|
pair.map {|r| sq * prod2 * r }
}
)
}
}
solutions.map {|pair| pair.sort } \
.uniq_by {|pair| pair[0] } \
.sort_by {|pair| pair[0] }
}
- Finds all solutions (a,b,c) such that: a^2 + b^2 + c^2 = n^2
func sum_of_three_squares(n) {
gather {
for k in (1 .. n//3) {
var t = sum_of_two_squares(n**2 - k**2) || next
take(t.map { [k, _...] }...)
}
}
}
say gather {
for n in (1..2200) {
sum_of_three_squares(n) || take(n)
}
}</lang>
- Output:
[1, 2, 4, 5, 8, 10, 16, 20, 32, 40, 64, 80, 128, 160, 256, 320, 512, 640, 1024, 1280, 2048]
Numbers d that cannot be expressed as a^2 + b^2 + c^2 = d^2, are numbers of the form 2^n or 5*2^n: <lang ruby>say gather {
for n in (1..2200) {
if ((n & (n-1) == 0) || (n%%5 && ((n/5) & (n/5 - 1) == 0))) {
take(n)
}
}
}</lang>
- Output:
[1, 2, 4, 5, 8, 10, 16, 20, 32, 40, 64, 80, 128, 160, 256, 320, 512, 640, 1024, 1280, 2048]
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