Pythagorean triples: Difference between revisions
(Add Seed7 example) |
(Updated D code) |
||
| Line 170: | Line 170: | ||
{{trans|C}} |
{{trans|C}} |
||
With the dmd compiler use -L/STACK:10000000 to increase stack size. |
With the dmd compiler use -L/STACK:10000000 to increase stack size. |
||
<lang d>import std.stdio |
<lang d>import std.stdio; |
||
// Should be ulong if going to or over 1 billion. |
// Should be ulong if going to or over 1 billion. |
||
| Line 178: | Line 178: | ||
void triples(in xint lim, |
void triples(in xint lim, |
||
in xint a=3, in xint b=4, in xint c=5) nothrow { |
in xint a=3, in xint b=4, in xint c=5) nothrow { |
||
immutable xint p = a + b + c; |
|||
[ 1, 2, 2, 2, 1, 2, 2, 2, 3], |
|||
[-1, 2, 2, -2, 1, 2, -2, 2, 3]]; |
|||
const xint p = a + b + c; |
|||
if (p > lim) |
if (p > lim) |
||
return; |
return; |
||
primitives++; |
primitives++; |
||
ntriples += lim / p; |
ntriples += lim / p; |
||
triples(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c); |
|||
triples(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c); |
|||
foreach (i; TypeTuple!(0, 1, 2)) |
|||
triples(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c); |
|||
triples(lim, |
|||
U[i][0] * a + U[i][1] * b + U[i][2] * c, |
|||
U[i][3] * a + U[i][4] * b + U[i][5] * c, |
|||
U[i][6] * a + U[i][7] * b + U[i][8] * c); |
|||
} |
} |
||
Revision as of 08:44, 20 August 2011
You are encouraged to solve this task according to the task description, using any language you may know.
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
Ada
Translation of efficient method from C, see the WP article. Compiles on gnat/gcc.
<lang Ada>with Ada.Text_IO;
procedure Pythagorean_Triples is
type Large_Natural is range 0 .. 2**63-1;
-- this is the maximum for gnat
procedure New_Triangle(A, B, C: Large_Natural;
Max_Perimeter: Large_Natural;
Total_Cnt, Primitive_Cnt: in out Large_Natural) is
Perimeter: constant Large_Natural := A + B + C;
begin
if Perimeter <= Max_Perimeter then
Primitive_Cnt := Primitive_Cnt + 1;
Total_Cnt := Total_Cnt + Max_Perimeter / Perimeter;
New_Triangle(A-2*B+2*C, 2*A-B+2*C, 2*A-2*B+3*C, Max_Perimeter, Total_Cnt, Primitive_Cnt);
New_Triangle(A+2*B+2*C, 2*A+B+2*C, 2*A+2*B+3*C, Max_Perimeter, Total_Cnt, Primitive_Cnt);
New_Triangle(2*B+2*C-A, B+2*C-2*A, 2*B+3*C-2*A, Max_Perimeter, Total_Cnt, Primitive_Cnt);
end if;
end New_Triangle;
T_Cnt, P_Cnt: Large_Natural;
begin
for I in 1 .. 9 loop
T_Cnt := 0;
P_Cnt := 0;
New_Triangle(3,4,5, 10**I, Total_Cnt => T_Cnt, Primitive_Cnt => P_Cnt);
Ada.Text_IO.Put_Line("Up to 10 **" & Integer'Image(I) & " :" &
Large_Natural'Image(T_Cnt) & " Triples," &
Large_Natural'Image(P_Cnt) & " Primitives");
end loop;
end Pythagorean_Triples;</lang>
Output:
Up to 10 ** 1 : 0 Triples, 0 Primitives Up to 10 ** 2 : 17 Triples, 7 Primitives Up to 10 ** 3 : 325 Triples, 70 Primitives Up to 10 ** 4 : 4858 Triples, 703 Primitives Up to 10 ** 5 : 64741 Triples, 7026 Primitives Up to 10 ** 6 : 808950 Triples, 70229 Primitives Up to 10 ** 7 : 9706567 Triples, 702309 Primitives Up to 10 ** 8 : 113236940 Triples, 7023027 Primitives Up to 10 ** 9 : 1294080089 Triples, 70230484 Primitives
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> Efficient method, generating primitive triples only as described in the same WP article:<lang C>#include <stdio.h>
- include <stdlib.h>
- include <stdint.h>
/* should be 64-bit integers if going over 1 billion */ typedef unsigned long xint;
- define FMT "%lu"
xint total, prim, max_peri; xint U[][9] = {{ 1, -2, 2, 2, -1, 2, 2, -2, 3}, { 1, 2, 2, 2, 1, 2, 2, 2, 3}, {-1, 2, 2, -2, 1, 2, -2, 2, 3}};
void new_tri(xint in[]) { int i; xint t[3], p = in[0] + in[1] + in[2];
if (p > max_peri) return;
prim ++;
/* for every primitive triangle, its multiples would be right-angled too; * count them up to the max perimeter */ total += max_peri / p;
/* recursively produce next tier by multiplying the matrices */ for (i = 0; i < 3; i++) { t[0] = U[i][0] * in[0] + U[i][1] * in[1] + U[i][2] * in[2]; t[1] = U[i][3] * in[0] + U[i][4] * in[1] + U[i][5] * in[2]; t[2] = U[i][6] * in[0] + U[i][7] * in[1] + U[i][8] * in[2]; new_tri(t); } }
int main() { xint seed[3] = {3, 4, 5};
for (max_peri = 10; max_peri <= 100000000; max_peri *= 10) { total = prim = 0; new_tri(seed);
printf( "Up to "FMT": "FMT" triples, "FMT" primitives.\n", max_peri, total, prim); } return 0; }</lang>Output<lang>Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives.</lang>
D
With the dmd compiler use -L/STACK:10000000 to increase stack size. <lang d>import std.stdio;
// Should be ulong if going to or over 1 billion. alias ulong xint; __gshared xint ntriples, primitives;
void triples(in xint lim,
in xint a=3, in xint b=4, in xint c=5) nothrow {
immutable xint p = a + b + c;
if (p > lim)
return;
primitives++;
ntriples += lim / p;
triples(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c);
triples(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c);
triples(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c);
}
void main() {
foreach (p; 1 .. 11) {
ntriples = primitives = 0;
triples(10L ^^ p);
writefln("Up to %11d: %11d triples, %9d primitives.",
10L ^^ p, ntriples, primitives);
}
}</lang>
Output:
Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives. Up to 1000000000: 1294080089 triples, 70230484 primitives. Up to 10000000000: 14557915466 triples, 702304875 primitives.
Runtime up to 10_000_000_000: about 72 seconds.
A shorter version: <lang d>import std.stdio, std.range, std.algorithm;
ulong[2] tri(in ulong lim, in ulong a=3, in ulong b=4, in ulong c=5) {
const ulong l = a + b + c; if (l > lim) return [0, 0]; ulong[2] r = [1, lim / l]; r[]+= tri(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c)[]; r[]+= tri(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c)[]; r[]+= tri(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c)[]; return r;
}
void main() {
foreach (peri; map!q{ 10 ^^ a }(iota(1, 9)))
writeln(peri, " ", tri(peri));
}</lang> Output:
10 [0, 0] 100 [7, 17] 1000 [70, 325] 10000 [703, 4858] 100000 [7026, 64741] 1000000 [70229, 808950] 10000000 [702309, 9706567] 100000000 [7023027, 113236940]
Euphoria
<lang euphoria>function tri(atom lim, sequence in)
sequence r
atom p
p = in[1] + in[2] + in[3]
if p > lim then
return {0, 0}
end if
r = {1, floor(lim / p)}
r += tri(lim, { in[1]-2*in[2]+2*in[3], 2*in[1]-in[2]+2*in[3], 2*in[1]-2*in[2]+3*in[3]})
r += tri(lim, { in[1]+2*in[2]+2*in[3], 2*in[1]+in[2]+2*in[3], 2*in[1]+2*in[2]+3*in[3]})
r += tri(lim, {-in[1]+2*in[2]+2*in[3], -2*in[1]+in[2]+2*in[3], -2*in[1]+2*in[2]+3*in[3]})
return r
end function
atom max_peri max_peri = 10 while max_peri <= 100000000 do
printf(1,"%d: ", max_peri)
? tri(max_peri, {3, 4, 5})
max_peri *= 10
end while</lang>
Output:
10: {0,0}
100: {7,17}
1000: {70,325}
10000: {703,4858}
100000: {7026,64741}
1000000: {70229,808950}
10000000: {702309,9706567}
100000000: {7023027,113236940}
Fortran
<lang fortran>module triples
implicit none
integer :: max_peri, prim, total
integer :: u(9,3) = reshape((/ 1, -2, 2, 2, -1, 2, 2, -2, 3, &
1, 2, 2, 2, 1, 2, 2, 2, 3, &
-1, 2, 2, -2, 1, 2, -2, 2, 3 /), &
(/ 9, 3 /))
contains
recursive subroutine new_tri(in)
integer, intent(in) :: in(:) integer :: i integer :: t(3), p
p = sum(in) if (p > max_peri) return
prim = prim + 1 total = total + max_peri / p do i = 1, 3 t(1) = sum(u(1:3, i) * in) t(2) = sum(u(4:6, i) * in) t(3) = sum(u(7:9, i) * in) call new_tri(t); end do
end subroutine new_tri end module triples
program Pythagorean
use triples implicit none
integer :: seed(3) = (/ 3, 4, 5 /) max_peri = 10 do total = 0 prim = 0 call new_tri(seed) write(*, "(a, i10, 2(i10, a))") "Up to", max_peri, total, " triples", prim, " primitives" if(max_peri == 100000000) exit max_peri = max_peri * 10 end do
end program Pythagorean</lang>
Output:
Up to 10 0 triples 0 primitives Up to 100 17 triples 7 primitives Up to 1000 325 triples 70 primitives Up to 10000 4858 triples 703 primitives Up to 100000 64741 triples 7026 primitives Up to 1000000 808950 triples 70229 primitives Up to 10000000 9706567 triples 702309 primitives Up to 100000000 113236940 triples 7023027 primitives
Icon and Unicon
This uses the elegant formula (#IV) from Formulas for generating Pythagorean triples
<lang Icon> link numbers 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>
numbers.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. I suspect there may be faster gcd algorithms that would speed this up.
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
/:~ ; <@(,.~ # {. 1:)@(*/~ 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.
That said, we do not actually have to generate all the triples, we just need to count them. Thus:
<lang j>trc=:3 :0
'm n'=. |:(#~ 1 = 2 | +/"1)(#~ >/"1) ,/ ,"0/~ }. i. <. %: y
<.y%+/"1 (#~ 1 = 2 +./@{. |:) (#~ y >: +/"1)m (-&*: ,. +:@* ,. +&*:) n
)</lang>
The result is a list of positive integers, one number for each primitive triple which fits within the limit, giving the number of triples which are multiples of that primitive triple whose perimeter is no greater than the limiting perimeter.
<lang> (#,+/)trc 1e8 7023027 113236940</lang>
But note that J's memory footprint reached 6.7GB during the computation, so to compute larger values the computation would have to be broken up into reasonable sized blocks.
Java
Brute force
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 = BigInteger.valueOf(100),
peri2 = periLimit.divide(BigInteger.valueOf(2)),
peri3 = periLimit.divide(BigInteger.valueOf(3));
for(BigInteger a = ONE; a.compareTo(peri3) < 0; a = a.add(ONE)){
BigInteger aa = a.multiply(a);
for(BigInteger b = a.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 = b.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).equals(ONE)){
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.
Brute force primitives with scaling
Pythagorean triples/Java/Brute force primitives
Parent/child
(with limited modification for saving a few BigInteger operations)
This can also go "forever" theoretically. Letting it go to another order of magnitude overflowed the stack on the computer this was tested on. This version also does not show the triples as it goes, it only counts them. <lang java5>import java.math.BigInteger;
public class Triples{
public static BigInteger LIMIT; public static final BigInteger TWO = BigInteger.valueOf(2); public static final BigInteger THREE = BigInteger.valueOf(3); public static final BigInteger FOUR = BigInteger.valueOf(4); public static final BigInteger FIVE = BigInteger.valueOf(5); public static long primCount = 0; public static long tripCount = 0;
//I don't know Japanese :p
public static void parChild(BigInteger a, BigInteger b, BigInteger c){
BigInteger perim = a.add(b).add(c);
if(perim.compareTo(LIMIT) > 0) return;
primCount++; tripCount += LIMIT.divide(perim).longValue();
BigInteger a2 = TWO.multiply(a), b2 = TWO.multiply(b), c2 = TWO.multiply(c),
c3 = THREE.multiply(c);
parChild(a.subtract(b2).add(c2),
a2.subtract(b).add(c2),
a2.subtract(b2).add(c3));
parChild(a.add(b2).add(c2),
a2.add(b).add(c2),
a2.add(b2).add(c3));
parChild(a.negate().add(b2).add(c2),
TWO.negate().multiply(a).add(b).add(c2),
TWO.negate().multiply(a).add(b2).add(c3));
}
public static void main(String[] args){
for(long i = 100; i <= 10000000; i*=10){
LIMIT = BigInteger.valueOf(i);
primCount = tripCount = 0;
parChild(THREE, FOUR, FIVE);
System.out.println(LIMIT + ": " + tripCount + " triples, " + primCount + " primitive.");
}
}
}</lang> Output:
100: 17 triples, 7 primitive. 1000: 325 triples, 70 primitive. 10000: 4858 triples, 703 primitive. 100000: 64741 triples, 7026 primitive. 1000000: 808950 triples, 70229 primitive. 10000000: 9706567 triples, 702309 primitive.
Modula-3
Note that this code only works on 64bit machines (where INTEGER is 64 bits). Modula-3 provides a LONGINT type, which is 64 bits on 32 bit systems, but there is a bug in the implementation apparently. <lang modula3>MODULE PyTriple64 EXPORTS Main;
IMPORT IO, Fmt;
VAR tcnt, pcnt, max, i: INTEGER;
PROCEDURE NewTriangle(a, b, c: INTEGER; VAR tcount, pcount: INTEGER) =
VAR perim := a + b + c;
BEGIN
IF perim <= max THEN
pcount := pcount + 1;
tcount := tcount + max DIV perim;
NewTriangle(a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c, tcount, pcount);
NewTriangle(a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c, tcount, pcount);
NewTriangle(2*b+2*c-a, b+2*c-2*a, 2*b+3*c-2*a, tcount, pcount);
END;
END NewTriangle;
BEGIN
i := 100;
REPEAT
max := i;
tcnt := 0;
pcnt := 0;
NewTriangle(3, 4, 5, tcnt, pcnt);
IO.Put(Fmt.Int(i) & ": " & Fmt.Int(tcnt) & " Triples, " &
Fmt.Int(pcnt) & " Primitives\n");
i := i * 10;
UNTIL i = 10000000;
END PyTriple64.</lang>
Output:
100: 17 Triples, 7 Primitives 1000: 325 Triples, 70 Primitives 10000: 4858 Triples, 703 Primitives 100000: 64741 Triples, 7026 Primitives 1000000: 808950 Triples, 70229 Primitives
Perl 6
First up is a fairly naive brute force implementation that is not really practical for large perimeter limits. 100 is no problem. 1000 is slow. 10000 is glacial.
Works with Rakudo 2011.06.
<lang perl6>my %triples; my $limit = 100;
for 3 .. $limit/2 -> $c {
for 1 .. $c -> $a {
my $b = ($c * $c - $a * $a).sqrt;
last if $c + $a + $b > $limit;
last if $a > $b;
if $b == $b.Int {
my $key = "$a $b $c";
%triples{$key} = ([gcd] $c, $a, $b.Int) > 1 ?? 0 !! 1;
say $key, %triples{$key} ?? ' - primitive' !! ;
}
}
}
say "There are {+%triples.keys} Pythagorean Triples with a perimeter <= $limit,"
~"\nof which {[+] %triples.values} are primitive.";</lang>
3 4 5 - primitive 6 8 10 5 12 13 - primitive 9 12 15 8 15 17 - primitive 12 16 20 7 24 25 - primitive 15 20 25 10 24 26 20 21 29 - primitive 18 24 30 16 30 34 21 28 35 12 35 37 - primitive 15 36 39 24 32 40 9 40 41 - primitive There are 17 Pythagorean Triples with a perimeter <= 100, of which 7 are primitive.
Here's a much faster version. Hint, "oyako" is Japanese for "parent/child". :-)
<lang perl6>sub triples($limit) {
my $primitive = 0;
my $civilized = 0;
sub oyako($a, $b, $c) {
my $perim = $a + $b + $c;
return if $perim > $limit;
++$primitive; $civilized += $limit div $perim;
oyako( $a - 2*$b + 2*$c, 2*$a - $b + 2*$c, 2*$a - 2*$b + 3*$c);
oyako( $a + 2*$b + 2*$c, 2*$a + $b + 2*$c, 2*$a + 2*$b + 3*$c);
oyako(-$a + 2*$b + 2*$c, -2*$a + $b + 2*$c, -2*$a + 2*$b + 3*$c);
}
oyako(3,4,5);
"$limit => ($primitive $civilized)";
}
for 10,100,1000 ... * -> $limit {
say triples $limit;
}</lang> Output:
10 => (0 0) 100 => (7 17) 1000 => (70 325) 10000 => (703 4858) 100000 => (7026 64741) 1000000 => (70229 808950) 10000000 => (702309 9706567) 100000000 => (7023027 113236940) 1000000000 => (70230484 1294080089) ^C
The geometric sequence of limits will continue on forever, so eventually when you get tired of waiting (about a billion on my computer), you can just stop it. Another efficiency trick of note: we avoid passing the limit in as a parameter to the inner helper routine, but instead supply the limit via the lexical scope. Likewise, the accumulators are referenced lexically, so only the triples themselves need to be passed downward, and nothing needs to be returned.
Here is an alternate version that avoids naming any scalars that can be handled by vector processing instead: <lang perl6>constant @coeff = [[+1, -2, +2], [+2, -1, +2], [+2, -2, +3]],
[[+1, +2, +2], [+2, +1, +2], [+2, +2, +3]],
[[-1, +2, +2], [-2, +1, +2], [-2, +2, +3]];
sub triples($limit) {
sub oyako(@trippy) {
my $perim = [+] @trippy;
return if $perim > $limit;
take (1 + ($limit div $perim)i);
for @coeff -> @nine {
oyako (map -> @three { [+] @three »*« @trippy }, @nine);
}
return;
}
my $complex = 0i + [+] gather oyako([3,4,5]);
"$limit => ({$complex.re, $complex.im})";
}
for 10,100,1000 ... * -> $limit {
say triples $limit;
}</lang> Using vectorized ops allows a bit more potential for parallelization, though this is probably not as big a win in this case, especially since we do a certain amount of multiplying by 1 that the scalar version doesn't need to do. Note the cute trick of adding complex numbers to add two numbers in parallel. The use of gather/take allows the summation to run in a different thread than the helper function, at least in theory...
In practice, this solution runs considerably slower than the previous one, due primarily to passing gather/take values up many levels of dynamic scope. Eventually this may be optimized. Also, this solution currently chews up gigabytes of memory, while the previous solution stays at a quarter gig or so. This likely indicates a memory leak somewhere in niecza.
PicoLisp
<lang PicoLisp>(for (Max 10 (>= 100000000 Max) (* Max 10))
(let (Total 0 Prim 0 In (3 4 5))
(recur (In)
(let P (apply + In)
(when (>= Max P)
(inc 'Prim)
(inc 'Total (/ Max P))
(for Row
(quote
(( 1 -2 2) ( 2 -1 2) ( 2 -2 3))
(( 1 2 2) ( 2 1 2) ( 2 2 3))
((-1 2 2) (-2 1 2) (-2 2 3)) )
(recurse
(mapcar '((U) (sum * U In)) Row) ) ) ) ) )
(prinl "Up to " Max ": " Total " triples, " Prim " primitives.") ) )</lang>
Output:
Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives.
Python
Two methods, the second of which is much faster <lang python>from fractions import gcd
def pt1(maxperimeter=100):
- Naive method
trips = []
for a in range(1, maxperimeter):
aa = a*a
for b in range(a, maxperimeter-a+1):
bb = b*b
for c in range(b, maxperimeter-b-a+1):
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 pytrip(trip=(3,4,5),perim=100, prim=1):
a0, b0, c0 = a, b, c = sorted(trip)
t, firstprim = set(), prim>0
while a + b + c <= perim:
t.add((a, b, c, firstprim>0))
a, b, c, firstprim = a+a0, b+b0, c+c0, False
#
t2 = set()
for a, b, c, firstprim in t:
a2, a5, b2, b5, c2, c3, c7 = a*2, a*5, b*2, b*5, c*2, c*3, c*7
if a5 - b5 + c7 <= perim:
t2 |= pytrip(( a - b2 + c2, a2 - b + c2, a2 - b2 + c3), perim, firstprim)
if a5 + b5 + c7 <= perim:
t2 |= pytrip(( a + b2 + c2, a2 + b + c2, a2 + b2 + c3), perim, firstprim)
if -a5 + b5 + c7 <= perim:
t2 |= pytrip((-a + b2 + c2, -a2 + b + c2, -a2 + b2 + c3), perim, firstprim)
return t | t2
def pt2(maxperimeter=100):
- Parent/child relationship method:
- http://en.wikipedia.org/wiki/Formulas_for_generating_Pythagorean_triples#XI.
trips = pytrip((3,4,5), maxperimeter, 1) return trips
def printit(maxperimeter=100, pt=pt1):
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 algo, mn, mx in ((pt1, 250, 2500), (pt2, 500, 20000)):
print(algo.__doc__)
for maxperimeter in range(mn, mx+1, mn):
printit(maxperimeter, algo)
</lang>
- Output
# Naive method
Up to a perimeter of 250 there are 56 triples, of which 18 are primitive
Up to a perimeter of 500 there are 137 triples, of which 35 are primitive
Up to a perimeter of 750 there are 227 triples, of which 52 are primitive
Up to a perimeter of 1000 there are 325 triples, of which 70 are primitive
Up to a perimeter of 1250 there are 425 triples, of which 88 are primitive
Up to a perimeter of 1500 there are 527 triples, of which 104 are primitive
Up to a perimeter of 1750 there are 637 triples, of which 123 are primitive
Up to a perimeter of 2000 there are 744 triples, of which 140 are primitive
Up to a perimeter of 2250 there are 858 triples, of which 156 are primitive
Up to a perimeter of 2500 there are 969 triples, of which 175 are primitive
# Parent/child relationship method:
# http://en.wikipedia.org/wiki/Formulas_for_generating_Pythagorean_triples#XI.
Up to a perimeter of 500 there are 137 triples, of which 35 are primitive
Up to a perimeter of 1000 there are 325 triples, of which 70 are primitive
Up to a perimeter of 1500 there are 527 triples, of which 104 are primitive
Up to a perimeter of 2000 there are 744 triples, of which 140 are primitive
Up to a perimeter of 2500 there are 969 triples, of which 175 are primitive
Up to a perimeter of 3000 there are 1204 triples, of which 211 are primitive
Up to a perimeter of 3500 there are 1443 triples, of which 245 are primitive
Up to a perimeter of 4000 there are 1687 triples, of which 280 are primitive
Up to a perimeter of 4500 there are 1931 triples, of which 314 are primitive
Up to a perimeter of 5000 there are 2184 triples, of which 349 are primitive
Up to a perimeter of 5500 there are 2442 triples, of which 385 are primitive
Up to a perimeter of 6000 there are 2701 triples, of which 422 are primitive
Up to a perimeter of 6500 there are 2963 triples, of which 457 are primitive
Up to a perimeter of 7000 there are 3224 triples, of which 492 are primitive
Up to a perimeter of 7500 there are 3491 triples, of which 527 are primitive
Up to a perimeter of 8000 there are 3763 triples, of which 560 are primitive
Up to a perimeter of 8500 there are 4029 triples, of which 597 are primitive
Up to a perimeter of 9000 there are 4304 triples, of which 631 are primitive
Up to a perimeter of 9500 there are 4577 triples, of which 667 are primitive
Up to a perimeter of 10000 there are 4858 triples, of which 703 are primitive
Up to a perimeter of 10500 there are 5138 triples, of which 736 are primitive
Up to a perimeter of 11000 there are 5414 triples, of which 770 are primitive
Up to a perimeter of 11500 there are 5699 triples, of which 804 are primitive
Up to a perimeter of 12000 there are 5980 triples, of which 839 are primitive
Up to a perimeter of 12500 there are 6263 triples, of which 877 are primitive
Up to a perimeter of 13000 there are 6559 triples, of which 913 are primitive
Up to a perimeter of 13500 there are 6843 triples, of which 949 are primitive
Up to a perimeter of 14000 there are 7129 triples, of which 983 are primitive
Up to a perimeter of 14500 there are 7420 triples, of which 1019 are primitive
Up to a perimeter of 15000 there are 7714 triples, of which 1055 are primitive
Up to a perimeter of 15500 there are 8004 triples, of which 1089 are primitive
Up to a perimeter of 16000 there are 8304 triples, of which 1127 are primitive
Up to a perimeter of 16500 there are 8595 triples, of which 1159 are primitive
Up to a perimeter of 17000 there are 8884 triples, of which 1192 are primitive
Up to a perimeter of 17500 there are 9189 triples, of which 1228 are primitive
Up to a perimeter of 18000 there are 9484 triples, of which 1264 are primitive
Up to a perimeter of 18500 there are 9791 triples, of which 1301 are primitive
Up to a perimeter of 19000 there are 10088 triples, of which 1336 are primitive
Up to a perimeter of 19500 there are 10388 triples, of which 1373 are primitive
Up to a perimeter of 20000 there are 10689 triples, of which 1408 are primitive
Barebone minimum for this task:<lang Python>from sys import setrecursionlimit setrecursionlimit(2000) # 2000 ought to be big enough for everybody
def triples(lim, a = 3, b = 4, c = 5): l = a + b + c if l > lim: return (0, 0) return reduce(lambda x, y: (x[0] + y[0], x[1] + y[1]), [ (1, lim / l), triples(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c), triples(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c), triples(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c) ])
for peri in [10 ** e for e in range(1, 8)]: print peri, triples(peri)</lang>Output:<lang>10 (0, 0) 100 (7, 17) 1000 (70, 325) 10000 (703, 4858) 100000 (7026, 64741) 1000000 (70229, 808950) 10000000 (702309, 9706567)</lang>
Seed7
The example below uses bigInteger numbers:
<lang seed7>$ include "seed7_05.s7i";
include "bigint.s7i";
var bigInteger: total is 0_; var bigInteger: prim is 0_; var bigInteger: max_peri is 10_;
const proc: new_tri (in bigInteger: a, in bigInteger: b, in bigInteger: c) is func
local
var bigInteger: p is 0_;
begin
p := a + b + c;
if p <= max_peri then
incr(prim);
total +:= max_peri div p;
new_tri( a - 2_*b + 2_*c, 2_*a - b + 2_*c, 2_*a - 2_*b + 3_*c);
new_tri( a + 2_*b + 2_*c, 2_*a + b + 2_*c, 2_*a + 2_*b + 3_*c);
new_tri(-a + 2_*b + 2_*c, -2_*a + b + 2_*c, -2_*a + 2_*b + 3_*c);
end if;
end func;
const proc: main is func
begin
while max_peri <= 100000000_ do
total := 0_;
prim := 0_;
new_tri(3_, 4_, 5_);
writeln("Up to " <& max_peri <& ": " <& total <& " triples, " <& prim <& " primitives.");
max_peri *:= 10_;
end while;
end func;</lang>
Output:
Up to 10: 0 triples, 0 primitives. Up to 100: 17 triples, 7 primitives. Up to 1000: 325 triples, 70 primitives. Up to 10000: 4858 triples, 703 primitives. Up to 100000: 64741 triples, 7026 primitives. Up to 1000000: 808950 triples, 70229 primitives. Up to 10000000: 9706567 triples, 702309 primitives. Up to 100000000: 113236940 triples, 7023027 primitives.
Tcl
Using the efficient method based off the Wikipedia article: <lang tcl>proc countPythagoreanTriples {limit} {
lappend q 3 4 5
set idx [set count [set prim 0]]
while {$idx < [llength $q]} {
set a [lindex $q $idx] set b [lindex $q [incr idx]] set c [lindex $q [incr idx]] incr idx if {$a + $b + $c <= $limit} { incr prim for {set i 1} {$i*$a+$i*$b+$i*$c <= $limit} {incr i} { incr count } lappend q \ [expr {$a + 2*($c-$b)}] [expr {2*($a+$c) - $b}] [expr {2*($a-$b) + 3*$c}] \ [expr {$a + 2*($b+$c)}] [expr {2*($a+$c) + $b}] [expr {2*($a+$b) + 3*$c}] \ [expr {2*($b+$c) - $a}] [expr {2*($c-$a) + $b}] [expr {2*($b-$a) + 3*$c}] }
} return [list $count $prim]
} for {set i 10} {$i <= 10000000} {set i [expr {$i*10}]} {
lassign [countPythagoreanTriples $i] count primitive puts "perimeter limit $i => $count triples, $primitive primitive"
}</lang> Output:
perimeter limit 10 => 0 triples, 0 primitive perimeter limit 100 => 17 triples, 7 primitive perimeter limit 1000 => 325 triples, 70 primitive perimeter limit 10000 => 4858 triples, 703 primitive perimeter limit 100000 => 64741 triples, 7026 primitive perimeter limit 1000000 => 808950 triples, 70229 primitive perimeter limit 10000000 => 9706567 triples, 702309 primitive