Pythagorean triples
A Pythagorean triple is defined as three positive integers where , and They are called primitive triples if are coprime, that is, if their pairwise greatest common divisors . Because of their relationship through the Pythagorean theorem, a, b, and c are coprime if a and b are coprime (). Each triple forms the length of the sides of a right triangle, whose perimeter is .
Task
The task is to determine how many Pythagorean triples there are with a perimeter no larger than 100 and the number of these that are primitive.
Extra credit: Deal with large values. Can your program handle a max perimeter of 1,000,000? What about 10,000,000? 100,000,000?
Note: the extra credit is not for you to demonstrate how fast your language is compared to others; you need a proper algorithm to solve them in a timely manner.
- Cf
C
Sample implemention; naive method, patentedly won't scale to larger numbers, despite the attempt to optimize it. Calculating up to 10000 is already a test of patience. <lang C>#include <stdio.h>
- include <stdlib.h>
typedef unsigned long long xint; typedef unsigned long ulong;
inline ulong gcd(ulong m, ulong n) { ulong t; while (n) { t = n; n = m % n; m = t; } return m; }
int main() { ulong a, b, c, pytha = 0, prim = 0, max_p = 100; xint aa, bb, cc;
for (a = 1; a <= max_p / 3; a++) { aa = (xint)a * a; printf("a = %lu\r", a); /* show that we are working */ fflush(stdout);
/* max_p/2: valid limit, because one side of triangle * must be less than the sum of the other two */ for (b = a + 1; b < max_p/2; b++) { bb = (xint)b * b; for (c = b + 1; c < max_p/2; c++) { cc = (xint)c * c; if (aa + bb < cc) break; if (a + b + c > max_p) break;
if (aa + bb == cc) { pytha++; if (gcd(a, b) == 1) prim++; } } } }
printf("Up to %lu, there are %lu triples, of which %lu are primitive\n", max_p, pytha, prim);
return 0; }</lang>output:<lang>Up to 100, there are 17 triples, of which 7 are primitive</lang>
Icon and Unicon
This uses the elegant formula (#IV) from Formulas for generating Pythagorean triples
<lang Icon> link math link printf
procedure main(A) # P-triples
plimit := (0 < integer(\A[1])) | 100 # get perimiter limit
nonprimitiveS := set() # record unique non-primitives triples
primitiveS := set() # record unique primitive triples
u := 0
while (g := (u +:= 1)^2) + 3 * u + 2 < plimit / 2 do {
every v := seq(1) do {
a := g + (i := 2*u*v)
b := (h := 2*v^2) + i
c := g + h + i
if (p := a + b + c) > plimit then break
insert( (gcd(u,v)=1 & u%2=1, primitiveS) | nonprimitiveS, memo(a,b,c))
every k := seq(2) do { # k is for larger non-primitives
if k*p > plimit then break
insert(nonprimitiveS,memo(a*k,b*k,c*k) )
}
}
}
printf("Under perimiter=%d: Pythagorean Triples=%d including primitives=%d\n",
plimit,*nonprimitiveS+*primitiveS,*primitiveS)
every put(gcol := [] , &collections) printf("Time=%d, Collections: total=%d string=%d block=%d",&time,gcol[1],gcol[3],gcol[4]) end
procedure memo(x[]) #: return a csv string of arguments in sorted order
every (s := "") ||:= !sort(x) do s ||:= ","
return s[1:-1]
end</lang>
math.icn provides gcd printf.icn provides printf
The output from some sample runs with BLKSIZE=500000000 and STRSIZE=50000000 is below. It starts getting very slow after 10M at about 9 minutes (times are shown in ms.
Output:
Under perimiter=10: Pythagorean Triples=0 including primitives=0 Time=3, Collections: total=0 string=0 block=0 Under perimiter=100: Pythagorean Triples=17 including primitives=7 Time=3, Collections: total=0 string=0 block=0 Under perimiter=1000: Pythagorean Triples=325 including primitives=70 Time=6, Collections: total=0 string=0 block=0 Under perimiter=10000: Pythagorean Triples=4858 including primitives=703 Time=57, Collections: total=0 string=0 block=0 Under perimiter=100000: Pythagorean Triples=64741 including primitives=7026 Time=738, Collections: total=0 string=0 block=0 Under perimiter=1000000: Pythagorean Triples=808950 including primitives=70229 Time=12454, Collections: total=0 string=0 block=0 Under perimiter=10000000: Pythagorean Triples=9706567 including primitives=702309 Time=560625, Collections: total=16 string=8 block=8
J
Brute force approach:
<lang j>pytr=:3 :0
r=.i.0 3
for_a.1+i.<.(y-1)%3 do.
b=.1+a+i.<.(y%2)-3*a%2
c=. a +&.*: b
keep=. (c=<.c)*.y>:a+b+c
if.1 e.keep do.
r=.r, a,.b ,.&(keep&#) c
end.
end.
(,.~ prim"1)r
)
prim=:1 = 2 +./@{. |:</lang>
Example use:
First column indicates whether the triple is primitive, and the remaining three columns are a, b and c.
<lang j> pytr 100 1 3 4 5 1 5 12 13 0 6 8 10 1 7 24 25 1 8 15 17 0 9 12 15 1 9 40 41 0 10 24 26 0 12 16 20 1 12 35 37 0 15 20 25 0 15 36 39 0 16 30 34 0 18 24 30 1 20 21 29 0 21 28 35 0 24 32 40
(# , [: {. +/) pytr 10
0 0
(# , [: {. +/) pytr 100
17 7
(# , [: {. +/) pytr 1000
325 70
(# , [: {. +/) pytr 10000
4858 703</lang>
pytr 10000 takes 4 seconds on this laptop, and time to complete grows with square of perimeter, so pytr 1e6 should take something like 11 hours using this algorithm on this machine.
A slightly smarter approach:
<lang j>trips=:3 :0
'm n'=. |:(#~ 1 = 2 | +/"1)(#~ >/"1),/,"0/~}.i.<.%:y
prim=. (#~ 1 = 2 +./@{. |:) (#~ y>:+/"1)m (-&*:,. +:@*,. +&*:) n
(,.~ e.&prim) ~./:~;<@(*/~ 1 + y i.@<.@% +/)"1 prim
)</lang>
usage for trips is the same as for pytr. Thus:
<lang j> (# , 1 {. +/) trips 10 0 0
(# , 1 {. +/) trips 100
17 7
(# , 1 {. +/) trips 1000
325 70
(# , 1 {. +/) trips 10000
4858 703
(# , 1 {. +/) trips 100000
64741 7026
(# , 1 {. +/) trips 1000000
808950 70229
(# , 1 {. +/) trips 10000000
9706567 702309</lang>
The last line took about 16 seconds.
Java
Theoretically, this can go "forever", but it takes a while, so only the minimum is shown. Luckily, BigInteger has a GCD method built in.
<lang java> import java.math.BigInteger; import static java.math.BigInteger.ONE;
public class PythTrip{
public static void main(String[] args){
long tripCount = 0, primCount = 0;
//change this to whatever perimeter limit you want;the RAM's the limit
BigInteger periLimit = new BigInteger("10000"),
peri2 = periLimit.divide(new BigInteger("2")),
peri3 = periLimit.divide(new BigInteger("3"));
for(BigInteger a = ONE; a.compareTo(peri3) < 0; a = a.add(ONE)){
BigInteger aa = a.multiply(a);
for(BigInteger b = new BigInteger(a.toString()).add(ONE);
b.compareTo(peri2) < 0; b = b.add(ONE)){
BigInteger bb = b.multiply(b);
BigInteger ab = a.add(b);
BigInteger aabb = aa.add(bb);
for(BigInteger c = new BigInteger(b.toString()).add(ONE);
c.compareTo(peri2) < 0; c = c.add(ONE)){
int compare = aabb.compareTo(c.multiply(c));
//if a+b+c > periLimit
if(ab.add(c).compareTo(periLimit) > 0){
break;
}
//if a^2 + b^2 != c^2
if(compare < 0){
break;
}else if (compare == 0){
tripCount++;
System.out.print(a + ", " + b + ", " + c);
//does binary GCD under the hood
if(a.gcd(b).compareTo(ONE) == 0){
System.out.print(" primitive");
primCount++;
}
System.out.println();
}
}
}
}
System.out.println("Up to a perimeter of " + periLimit + ", there are "
+ tripCount + " triples, of which " + primCount + " are primitive.");
}
}</lang> Output:
3, 4, 5 primitive 5, 12, 13 primitive 6, 8, 10 7, 24, 25 primitive 8, 15, 17 primitive 9, 12, 15 9, 40, 41 primitive 10, 24, 26 12, 16, 20 12, 35, 37 primitive 15, 20, 25 15, 36, 39 16, 30, 34 18, 24, 30 20, 21, 29 primitive 21, 28, 35 24, 32, 40 Up to a perimeter of 100, there are 17 triples, of which 7 are primitive.
Python
I could not wait for the answer to triples with a perimeter of 10000, although it should calculate it. <lang python>>>> from fractions import gcd >>> def pt(maxperimeter=100): trips = [] for a in range(1, maxperimeter): aa = a*a for b in range(a, maxperimeter-a): bb = b*b for c in range(b, maxperimeter-b-a): cc = c*c if a+b+c > maxperimeter or cc > aa + bb: break if aa + bb == cc: trips.append((a,b,c, gcd(a, b) == 1)) return trips
>>> def printit(maxperimeter=100): trips = pt(maxperimeter) print("Up to a perimeter of %i there are %i triples, of which %i are primitive" % (maxperimeter, len(trips), len([prim for a,b,c,prim in trips if prim])))
>>> for i in range(4): printit(10**i)
Up to a perimeter of 1 there are 0 triples, of which 0 are primitive Up to a perimeter of 10 there are 0 triples, of which 0 are primitive Up to a perimeter of 100 there are 17 triples, of which 7 are primitive Up to a perimeter of 1000 there are 324 triples, of which 70 are primitive</lang>