Pythagorean triples: Difference between revisions
m (→{{header|REXX}}: optimized routine for finding a Pythagorean triple.) |
(→{{header|Perl 6}}: update) |
||
| Line 2,237: | Line 2,237: | ||
=={{header|Perl 6}}== |
=={{header|Perl 6}}== |
||
| ⚫ | |||
Here is a straight-forward, naive brute force implementation: |
Here is a straight-forward, naive brute force implementation: |
||
<lang perl6>constant limit = |
<lang perl6>constant limit = 100; |
||
for [X] [^limit] xx 3 -> \a, \b, \c { |
for [X] [^limit] xx 3 -> (\a, \b, \c) { |
||
say [a, b, c] if a < b < c and a**2 + b**2 == c**2 |
say [a, b, c] if a < b < c and a**2 + b**2 == c**2 |
||
}</lang> |
}</lang> |
||
Here is a slightly less naive brute force implementation that is not really practical for large perimeter limits. |
Here is a slightly less naive brute force implementation that is not really practical for large perimeter limits. It's pretty zippy up to about 10000 though. |
||
<lang perl6>my %triples; |
<lang perl6>my %triples; |
||
my $limit = |
my $limit = 10000; |
||
for 3 .. $limit/2 -> $c { |
for 3 .. $limit/2 -> $c { |
||
| Line 2,265: | Line 2,265: | ||
say "There are {+%triples.keys} Pythagorean Triples with a perimeter <= $limit," |
say "There are {+%triples.keys} Pythagorean Triples with a perimeter <= $limit," |
||
~"\nof which {[+] %triples.values} are primitive.";</lang> |
~"\nof which {[+] %triples.values} are primitive.";</lang> |
||
{{out}} |
|||
<pre>3 4 5 - primitive |
<pre>3 4 5 - primitive |
||
6 8 10 |
6 8 10 |
||
| Line 2,276: | Line 2,276: | ||
10 24 26 |
10 24 26 |
||
20 21 29 - primitive |
20 21 29 - primitive |
||
... |
|||
18 24 30 |
|||
196 4800 4804 |
|||
16 30 34 |
|||
310 4800 4810 |
|||
21 28 35 |
|||
392 4794 4810 |
|||
| ⚫ | |||
171 4872 4875 |
|||
15 36 39 |
|||
| ⚫ | |||
24 32 40 |
|||
140 4899 4901 - primitive |
|||
There are |
There are 4858 Pythagorean Triples with a perimeter <= 10000, |
||
of which |
of which 703 are primitive.</pre> |
||
</pre> |
|||
Here's a much faster version. Hint, "oyako" is Japanese for "parent/child". <tt>:-)</tt> |
Here's a much faster version. Hint, "oyako" is Japanese for "parent/child". <tt>:-)</tt> |
||
| ⚫ | |||
<lang perl6>sub triples($limit) { |
<lang perl6>sub triples($limit) { |
||
my $primitive = 0; |
my $primitive = 0; |
||
| Line 2,349: | Line 2,347: | ||
The use of <tt>gather</tt>/<tt>take</tt> allows the summation to run in a different thread than the helper function, at least in theory... |
The use of <tt>gather</tt>/<tt>take</tt> 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 <tt>gather</tt>/<tt>take</tt> values up many levels of dynamic scope. Eventually this may be optimized |
In practice, this solution runs considerably slower than the previous one, due primarily to passing <tt>gather</tt>/<tt>take</tt> values up many levels of dynamic scope. Eventually this may be optimized. |
||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
Revision as of 07:47, 21 October 2015
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
AutoHotkey
<lang autohotkey>#NoEnv SetBatchLines, -1
- SingleInstance, Force
- Greatest common divisor, from http://rosettacode.org/wiki/Greatest_common_divisor#AutoHotkey
gcd(a,b) { Return b=0 ? Abs(a) : Gcd(b,mod(a,b)) }
count_triples(max) { primitives := 0, triples := 0, m := 2 while m <= (max / 2)**0.5 { n := mod(m, 2) + 1 ,p := 2*m*(m + n) , delta := 4*m while n < m and p <= max gcd(m, n) = 1 ? (primitives++ , triples += max // p) : "" , n += 2 , p += delta m++ } Return primitives " primitives out of " triples " triples" }
Loop, 8 Msgbox % 10**A_Index ": " count_triples(10**A_Index)</lang>
- Output:
10: 0 primitives out of 0 triples 100: 7 primitives out of 17 triples 1000: 70 primitives out of 325 triples 10000: 703 primitives out of 4858 triples 100000: 7026 primitives out of 64741 triples 1000000: 70229 primitives out of 808950 triples 10000000: 702309 primitives out of 9706567 triples 100000000: 7023027 primitives out of 113236940 triples
BBC BASIC
The built-in array arithmetic is very well suited to this task! <lang bbcbasic> DIM U0%(2,2), U1%(2,2), U2%(2,2), seed%(2)
U0%() = 1, -2, 2, 2, -1, 2, 2, -2, 3
U1%() = 1, 2, 2, 2, 1, 2, 2, 2, 3
U2%() = -1, 2, 2, -2, 1, 2, -2, 2, 3
seed%() = 3, 4, 5
FOR power% = 1 TO 7
all% = 0 : prim% = 0
PROCtri(seed%(), 10^power%, all%, prim%)
PRINT "Up to 10^"; power%, ": " all% " triples" prim% " primitives"
NEXT
END
DEF PROCtri(i%(), mp%, RETURN all%, RETURN prim%)
LOCAL t%() : DIM t%(2)
IF SUM(i%()) > mp% ENDPROC
prim% += 1
all% += mp% DIV SUM(i%())
t%() = U0%() . i%()
PROCtri(t%(), mp%, all%, prim%)
t%() = U1%() . i%()
PROCtri(t%(), mp%, all%, prim%)
t%() = U2%() . i%()
PROCtri(t%(), mp%, all%, prim%)
ENDPROC</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
Bracmat
<lang bracmat>(pythagoreanTriples=
total prim max-peri U
. (.(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))
: ?U
& ( new-tri
= i t p Urows Urow Ucols
, a b c loop A B C
. !arg:(,?a,?b,?c)
& !a+!b+!c:~>!max-peri:?p
& 1+!prim:?prim
& div$(!max-peri.!p)+!total:?total
& !U:?Urows
& ( loop
= !Urows:(.?Urow) ?Urows
& !Urow:?Ucols
& :?t
& whl
' ( !Ucols:(?A,?B,?C) ?Ucols
& (!t,!a*!A+!b*!B+!c*!C):?t
)
& new-tri$!t
& !loop
)
& !loop
|
)
& ( Main
= seed
. (,3,4,5):?seed
& 10:?max-peri
& whl
' ( 0:?total:?prim
& new-tri$!seed
& out
$ ( str
$ ( "Up to "
!max-peri
": "
!total
" triples, "
!prim
" primitives."
)
)
& !max-peri*10:~>10000000:?max-peri
)
)
& Main$
);
pythagoreanTriples$; </lang>
Output (under Linux):
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.
Under Windows XP Command prompt the last result is unattainable due to stack overflow.
With very few changes we can get rid of the stack exhausting recursion. Instead of calling new-tri recursively, be push the triples to test onto a stack and return to the Main function. In the innermost loop we pop a triple from the stack and call new-tri. The memory overhead is only a few megabytes for a max perimeter of 100,000,000. On my Windows XP box the whole computation takes at least 15 minutes! Given enough time (and memory), the program can compute results for larger perimeters.
<lang bracmat>(pythagoreanTriples=
total prim max-peri U stack
. (.(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))
: ?U
& ( new-tri
= i t p Urows Urow Ucols Ucol
, a b c loop A B C
. !arg:(,?a,?b,?c)
& !a+!b+!c:~>!max-peri:?p
& 1+!prim:?prim
& div$(!max-peri.!p)+!total:?total
& !U:?Urows
& ( loop
= !Urows:(.?Urow) ?Urows
& !Urow:?Ucols
& :?t
& whl
' ( !Ucols:(?A,?B,?C) ?Ucols
& (!t,!a*!A+!b*!B+!c*!C):?t
)
& !t !stack:?stack
& !loop
)
& !loop
|
)
& ( Main
= seed
. 10:?max-peri
& whl
' ( 0:?total:?prim
& (,3,4,5):?stack
& whl
' (!stack:%?seed ?stack&new-tri$!seed)
& out
$ ( str
$ ( "Up to "
!max-peri
": "
!total
" triples, "
!prim
" primitives."
)
)
& !max-peri*10:~>100000000:?max-peri
)
)
& Main$
);
pythagoreanTriples$;</lang>
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>
Same as above, but with loop unwound and third recursion eliminated: <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;
void new_tri(xint in[]) {
int i; xint t[3], p; xint x = in[0], y = in[1], z = in[2];
recur: p = x + y + z;
if (p > max_peri) return;
prim ++; total += max_peri / p;
t[0] = x - 2 * y + 2 * z; t[1] = 2 * x - y + 2 * z; t[2] = t[1] - y + z; new_tri(t);
t[0] += 4 * y; t[1] += 2 * y; t[2] += 4 * y; new_tri(t);
z = t[2] - 4 * x; y = t[1] - 4 * x; x = t[0] - 2 * x; goto recur;
}
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>
C#
Based on Ada example, which is a translation of efficient method from C, see the WP article.
<lang C sharp>using System;
namespace RosettaCode.CSharp {
class Program
{
static void Count_New_Triangle(ulong A, ulong B, ulong C, ulong Max_Perimeter, ref ulong Total_Cnt, ref ulong Primitive_Cnt)
{
ulong Perimeter = A + B + C;
if (Perimeter <= Max_Perimeter)
{
Primitive_Cnt = Primitive_Cnt + 1;
Total_Cnt = Total_Cnt + Max_Perimeter / Perimeter;
Count_New_Triangle(A + 2 * C - 2 * B, 2 * A + 2 * C - B, 2 * A + 3 * C - 2 * B, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt);
Count_New_Triangle(A + 2 * B + 2 * C, 2 * A + B + 2 * C, 2 * A + 2 * B + 3 * C, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt);
Count_New_Triangle(2 * B + 2 * C - A, B + 2 * C - 2 * A, 2 * B + 3 * C - 2 * A, Max_Perimeter, ref Total_Cnt, ref Primitive_Cnt);
}
}
static void Count_Pythagorean_Triples()
{
ulong T_Cnt, P_Cnt;
for (int I = 1; I <= 8; I++)
{
T_Cnt = 0;
P_Cnt = 0;
ulong ExponentNumberValue = (ulong)Math.Pow(10, I);
Count_New_Triangle(3, 4, 5, ExponentNumberValue, ref T_Cnt, ref P_Cnt);
Console.WriteLine("Perimeter up to 10E" + I + " : " + T_Cnt + " Triples, " + P_Cnt + " Primitives");
}
}
static void Main(string[] args)
{
Count_Pythagorean_Triples();
}
}
}</lang>
Output:
Perimeter up to 10E1 : 0 Triples, 0 Primitives Perimeter up to 10E2 : 17 Triples, 7 Primitives Perimeter up to 10E3 : 325 Triples, 70 Primitives Perimeter up to 10E4 : 4858 Triples, 703 Primitives Perimeter up to 10E5 : 64741 Triples, 7026 Primitives Perimeter up to 10E6 : 808950 Triples, 70229 Primitives Perimeter up to 10E7 : 9706567 Triples, 702309 Primitives Perimeter up to 10E8 : 113236940 Triples, 7023027 Primitives
Clojure
This version is based on Euclid's formula: for each pair (m,n) such that m>n>0, m and n coprime and of opposite polarity (even/odd), there is a primitive Pythagorean triple. It can be proven that the converse is true as well. <lang clojure>(defn gcd [a b] (if (zero? b) a (recur b (mod a b))))
(defn pyth [peri]
(for [m (range 2 (Math/sqrt (/ peri 2)))
n (range (inc (mod m 2)) m 2) ; n<m, opposite polarity
:let [p (* 2 m (+ m n))] ; = a+b+c for this (m,n)
:while (<= p peri)
:when (= 1 (gcd m n))
:let [m2 (* m m), n2 (* n n),
[a b] (sort [(- m2 n2) (* 2 m n)]), c (+ m2 n2)]
k (range 1 (inc (quot peri p)))]
[(= k 1) (* k a) (* k b) (* k c)]))
(defn rcount [ts] ; (->> peri pyth rcount) produces [total, primitive] counts
(reduce (fn [[total prims] t] [(inc total), (if (first t) (inc prims) prims)]) [0 0] ts))</lang>
To handle really large perimeters, we can dispense with actually generating the triples and just calculate the counts: <lang clojure>(defn pyth-count [peri]
(reduce (fn [[total prims] k] [(+ total k), (inc prims)]) [0 0]
(for [m (range 2 (Math/sqrt (/ peri 2)))
n (range (inc (mod m 2)) m 2) ; n<m, opposite polarity
:let [p (* 2 m (+ m n))] ; = a+b+c for this (m,n)
:while (<= p peri)
:when (= 1 (gcd m n))]
(quot peri p))))</lang>
CoffeeScript
This algorithm scales linearly with the max perimeter. It uses two loops that are capped by the square root of the half-perimeter to examine/count provisional values of m and n, where m and n generate a, b, c, and p using simple number theory.
<lang coffeescript> gcd = (x, y) ->
return x if y == 0 gcd(y, x % y)
- m,n generate primitive Pythag triples
- preconditions:
- m, n are integers of different parity
- m > n
- gcd(m,n) == 1 (coprime)
- m, n generate: [m*m - n*n, 2*m*n, m*m + n*n]
- perimeter is 2*m*m + 2*m*n = 2 * m * (m+n)
count_triples = (max_perim) ->
num_primitives = 0
num_triples = 0
m = 2
upper_limit = Math.sqrt max_perim / 2
while m <= upper_limit
n = m % 2 + 1
p = 2*m*m + 2*m*n
delta = 4*m
while n < m and p <= max_perim
if gcd(m, n) == 1
num_primitives += 1
num_triples += Math.floor max_perim / p
n += 2
p += delta
m += 1
console.log num_primitives, num_triples
max_perim = Math.pow 10, 9 # takes under a minute count_triples(max_perim) </lang> output
time coffee pythag_triples.coffee 70230484 1294080089 real 0m45.989s
Common Lisp
<lang lisp>(defun mmul (a b)
(loop for x in a collect
(loop for y in x
for z in b sum (* y z))))
(defun count-tri (lim)
(let ((prim 0) (cnt 0))
(labels ((count1 (tr)
(let ((peri (reduce #'+ tr)))
(when (<= peri lim)
(incf prim)
(incf cnt (truncate lim peri))
(count1 (mmul '(( 1 -2 2) ( 2 -1 2) ( 2 -2 3)) tr))
(count1 (mmul '(( 1 2 2) ( 2 1 2) ( 2 2 3)) tr))
(count1 (mmul '((-1 2 2) (-2 1 2) (-2 2 3)) tr))))))
(count1 '(3 4 5))
(format t "~a: ~a prim, ~a all~%" lim prim cnt))))
(loop for p from 2 do (count-tri (expt 10 p)))</lang>output<lang>100: 7 prim, 17 all 1000: 70 prim, 325 all 10000: 703 prim, 4858 all 100000: 7026 prim, 64741 all 1000000: 70229 prim, 808950 all 10000000: 702309 prim, 9706567 all ...</lang>
D
Lazy Functional Version
With hints from the Haskell solution. <lang d>void main() @safe {
import std.stdio, std.range, std.algorithm, std.typecons, std.numeric;
enum triples = (in uint n) pure nothrow @safe /*@nogc*/ =>
iota(1, n + 1)
.map!(z => iota(1, z + 1)
.map!(x => iota(x, z + 1).map!(y => tuple(x, y, z))))
.joiner.joiner
.filter!(t => t[0] ^^ 2 + t[1] ^^ 2 == t[2] ^^ 2 && t[].only.sum <= n)
.map!(t => tuple(t[0 .. 2].gcd == 1, t[]));
auto xs = triples(100);
writeln("Up to 100 there are ", xs.count, " triples, ",
xs.filter!q{ a[0] }.count, " are primitive.");
}</lang>
- Output:
Up to 100 there are 17 triples, 7 are primitive.
Shorter Version
<lang d>ulong[2] tri(ulong lim, ulong a=3, ulong b=4, ulong c=5) pure nothrow @safe @nogc {
immutable l = a + b + c;
if (l > lim)
return [0, 0];
typeof(return) 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() /*@safe*/ {
import std.stdio;
foreach (immutable p; 1 .. 9)
writeln(10 ^^ p, ' ', tri(10 ^^ p));
}</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]
Run-time (32 bit system): about 0.80 seconds with ldc2.
Short SIMD Version
With LDC compiler this is a little faster than the precedent version (remove @nogc to compile it with the current version of LDC compiler). <lang d>import std.stdio, core.simd;
ulong2 tri(in ulong lim, in ulong a=3, in ulong b=4, in ulong c=5) pure nothrow @safe @nogc {
immutable l = a + b + c;
if (l > lim)
return [0, 0];
typeof(return) 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() /*@safe*/ {
foreach (immutable p; 1 .. 9)
writeln(10 ^^ p, ' ', tri(10 ^^ p).array);
}</lang> The output is the same. Run-time (32 bit system): about 0.67 seconds with ldc2.
Faster Version
<lang d>import std.stdio;
alias Xuint = uint; // ulong if going over 1 billion.
__gshared Xuint nTriples, nPrimitives, limit;
void countTriples(Xuint x, Xuint y, Xuint z) nothrow @nogc {
while (true) {
immutable p = x + y + z;
if (p > limit)
return;
nPrimitives++;
nTriples += limit / p;
auto t0 = x - 2 * y + 2 * z;
auto t1 = 2 * x - y + 2 * z;
auto t2 = t1 - y + z;
countTriples(t0, t1, t2);
t0 += 4 * y;
t1 += 2 * y;
t2 += 4 * y;
countTriples(t0, t1, t2);
z = t2 - 4 * x;
y = t1 - 4 * x;
x = t0 - 2 * x;
}
}
void main() {
foreach (immutable p; 1 .. 9) {
limit = Xuint(10) ^^ p;
nTriples = nPrimitives = 0;
countTriples(3, 4, 5);
writefln("Up to %11d: %11d triples, %9d primitives.",
limit, nTriples, nPrimitives);
}
}</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.
Run-time: about 0.27 seconds with ldc2.
Using the power p up to 11, using ulong for xuint, and compiling with the dmd -L/STACK:10000000 switch to increase the stack size to about 10MB:
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.
Total run-time up to 10_000_000_000: about 63 seconds.
Waiting less than half an hour:
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. Up to 100000000000: 161750315680 triples, 7023049293 primitives.
Eiffel
<lang Eiffel> class APPLICATION
create make
feature
make local perimeter: INTEGER do perimeter := 100 from until perimeter > 1000000 loop total := 0 primitive_triples := 0 count_pythagorean_triples (3, 4, 5, perimeter) io.put_string ("There are " + total.out + " triples, below " + perimeter.out + ". Of which " + primitive_triples.out + " are primitives.%N") perimeter := perimeter * 10 end end
count_pythagorean_triples (a, b, c, perimeter: INTEGER) -- Total count of pythagorean triples and total count of primitve triples below perimeter. local p: INTEGER do p := a + b + c if p <= perimeter then primitive_triples := primitive_triples + 1 total := total + perimeter // p count_pythagorean_triples (a + 2 * (- b + c), 2 * (a + c) - b, 2 * (a - b + c) + c, perimeter) count_pythagorean_triples (a + 2 * (b + c), 2 * (a + c) + b, 2 * (a + b + c) + c, perimeter) count_pythagorean_triples (- a + 2 * (b + c), 2 * (- a + c) + b, 2 * (- a + b + c) + c, perimeter) end end
feature {NONE}
primitive_triples: INTEGER
total: INTEGER
end </lang>
- Output:
There are 17 triples, below 100. Of which 7 are primitives. There are 325 triples, below 1000. Of which 70 are primitives. There are 4858 triples, below 10000. Of which 703 are primitives. There are 64741 triples, below 100000. Of which 7026 are primitives. There are 808950 triples, below 1000000. Of which 70229 are primitives.
Elixir
<lang elixir>defmodule RC do
def count_triples(limit), do: count_triples(limit,3,4,5)
defp count_triples(limit, a, b, c) when limit<(a+b+c), do: {0,0}
defp count_triples(limit, a, b, c) do
{p1, t1} = count_triples(limit, a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c)
{p2, t2} = count_triples(limit, a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c)
{p3, t3} = count_triples(limit,-a+2*b+2*c,-2*a+b+2*c,-2*a+2*b+3*c)
{1+p1+p2+p3, div(limit, a+b+c)+t1+t2+t3}
end
end
list = for n <- 1..8, do: Enum.reduce(1..n, 1, fn(_,acc)->10*acc end) Enum.each(list, fn n -> IO.inspect {n, RC.count_triples(n)} end)</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}}
Erlang
<lang Erlang>%% %% Pythagorian triples in Erlang, J.W. Luiten %% -module(triples). -export([main/1]).
%% Transformations t1, t2 and t3 to generate new triples t1(A, B, C) ->
{A-2*B+2*C, 2*A-B+2*C, 2*A-2*B+3*C}.
t2(A, B, C) ->
{A+2*B+2*C, 2*A+B+2*C, 2*A+2*B+3*C}.
t3(A, B, C) ->
{2*B+2*C-A, B+2*C-2*A, 2*B+3*C-2*A}.
%% Generation of triples count_triples(A, B, C, Tot_acc, Cnt_acc, Max_perimeter) when (A+B+C) =< Max_perimeter ->
Tot1 = Tot_acc + Max_perimeter div (A+B+C),
{A1, B1, C1} = t1(A, B, C),
{Tot2, Cnt2} = count_triples(A1, B1, C1, Tot1, Cnt_acc+1, Max_perimeter),
{A2, B2, C2} = t2(A, B, C),
{Tot3, Cnt3} = count_triples(A2, B2, C2, Tot2, Cnt2, Max_perimeter),
{A3, B3, C3} = t3(A, B, C),
{Tot4, Cnt4} = count_triples(A3, B3, C3, Tot3, Cnt3, Max_perimeter),
{Tot4, Cnt4};
count_triples(_A, _B, _C, Tot_acc, Cnt_acc, _Max_perimeter) ->
{Tot_acc, Cnt_acc}.
count_triples(A, B, C, Pow) ->
Max = trunc(math:pow(10, Pow)),
{Tot, Prim} = count_triples(A, B, C, 0, 0, Max),
{Pow, Tot, Prim}.
count_triples(Pow) ->
count_triples(3, 4, 5, Pow).
%% Display a single result. display_result({Pow, Tot, Prim}) ->
io:format("Up to 10 ** ~w : ~w triples, ~w primitives~n", [Pow, Tot, Prim]).
main(Max) ->
L = lists:seq(1, Max), Answer = lists:map(fun(X) -> count_triples(X) end, L), lists:foreach(fun(Result) -> display_result(Result) end, Answer).</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 Up to 10 ** 10 : 14557915466 triples, 702304875 primitives Up to 10 ** 11 : 161750315680 triples, 7023049293 primitives
ERRE
<lang ERRE>PROGRAM PIT
BEGIN
PRINT(CHR$(12);) !CLS PRINT(TIME$)
FOR POWER=1 TO 7 DO PLIMIT=10#^POWER UPPERBOUND=INT(1+PLIMIT^0.5) PRIMITIVES=0 TRIPLES=0 EXTRAS=0 ! will count the in-range multiples of any primitive
FOR M=2 TO UPPERBOUND DO
FOR N=1+(M MOD 2=1) TO M-1 STEP 2 DO
TERM1=2*M*N
TERM2=M*M-N*N
TERM3=M*M+N*N
PERIMETER=TERM1+TERM2+TERM3
IF PERIMETER<=PLIMIT THEN TRIPLES=TRIPLES+1
A=TERM1
B=TERM2
REPEAT
R=A-B*INT(A/B)
A=B
B=R
UNTIL R<=0
! we've found a primitive triple if a = 1, since hcf =1.
! and it is inside perimeter range. Save it in an array
IF (A=1) AND (PERIMETER<=PLIMIT) THEN
PRIMITIVES=PRIMITIVES+1
!-----------------------------------------------
!swap so in increasing order of side length
!-----------------------------------------------
IF TERM1>TERM2 THEN SWAP(TERM1,TERM2)
!-----------------------------------------------
!we have the primitive & removed any multiples.
!Now calculate ALL the multiples in range.
!-----------------------------------------------
NEX=INT(PLIMIT/PERIMETER)
EXTRAS=EXTRAS+NEX
END IF
!scan
END FOR
END FOR
PRINT("Primit. with perimeter <=";10#^power;"is";primitives;"&";extras;"non-prim.triples.")
PRINT(TIME$)
END FOR
PRINT PRINT("** End **")
END PROGRAM</lang>
- Output:
16:08:39 Primit. with perimeter <= 10 is 0 & 0 non-prim.triples. 16:08:39 Primit. with perimeter <= 100 is 7 & 17 non-prim.triples. 16:08:39 Primit. with perimeter <= 1000 is 70 & 325 non-prim.triples. 16:08:39 Primit. with perimeter <= 10000 is 703 & 4858 non-prim.triples. 16:08:39 Primit. with perimeter <= 100000 is 7026 & 64741 non-prim.triples. 16:08:41 Primit. with perimeter <= 1000000 is 70229 & 808950 non-prim.triples. 16:09:07 Primit. with perimeter <= 10000000 is 702309 & 9706567 non-prim.triples. 16:13:10 ** End **
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}
F#
<lang fsharp>let isqrt n =
let rec iter t =
let d = n - t*t
if (0 <= d) && (d < t+t+1) // t*t <= n < (t+1)*(t+1)
then t else iter ((t+(n/t))/2)
iter 1
let rec gcd a b =
let t = a % b if t = 0 then b else gcd b t
let coprime a b = gcd a b = 1
let num_to ms =
let mutable ctr = 0
let mutable prim_ctr = 0
let max_m = isqrt (ms/2)
for m = 2 to max_m do
for j = 0 to (m/2) - 1 do
let n = m-(2*j+1)
if coprime m n then
let s = 2*m*(m+n)
if s <= ms then
ctr <- ctr + (ms/s)
prim_ctr <- prim_ctr + 1
(ctr, prim_ctr)
let show i =
let s, p = num_to i in printfn "For perimeters up to %d there are %d total and %d primitive" i s p;;
List.iter show [ 100; 1000; 10000; 100000; 1000000; 10000000; 100000000 ]</lang>
- Output:
For perimeters up to 100 there are 17 total and 7 primitive For perimeters up to 1000 there are 325 total and 70 primitive For perimeters up to 10000 there are 4858 total and 703 primitive For perimeters up to 100000 there are 64741 total and 7026 primitive For perimeters up to 1000000 there are 808950 total and 70229 primitive For perimeters up to 10000000 there are 9706567 total and 702309 primitive For perimeters up to 100000000 there are 113236940 total and 7023027 primitive
Factor
Pretty slow (100 times slower than C)...
<lang factor>USING: accessors arrays formatting kernel literals math math.functions math.matrices math.ranges sequences ; IN: rosettacode.pyth
CONSTANT: T1 {
{ 1 2 2 }
{ -2 -1 -2 }
{ 2 2 3 }
} CONSTANT: T2 {
{ 1 2 2 }
{ 2 1 2 }
{ 2 2 3 }
} CONSTANT: T3 {
{ -1 -2 -2 }
{ 2 1 2 }
{ 2 2 3 }
}
CONSTANT: base { 3 4 5 }
TUPLE: triplets-count primitives total ;
- <0-triplets-count> ( -- a ) 0 0 \ triplets-count boa ;
- next-triplet ( triplet T -- triplet' ) [ 1array ] [ m. ] bi* first ;
- candidates-triplets ( seed -- candidates )
${ T1 T2 T3 } [ next-triplet ] with map ;
- add-triplets ( current-triples limit triplet -- stop )
sum 2dup > [ /i [ + ] curry change-total [ 1 + ] change-primitives drop t ] [ 3drop f ] if ;
- all-triplets ( current-triples limit seed -- triplets )
3dup add-triplets [ candidates-triplets [ all-triplets ] with swapd reduce ] [ 2drop ] if ;
- count-triplets ( limit -- count )
<0-triplets-count> swap base all-triplets ;
- pprint-triplet-count ( limit count -- )
[ total>> ] [ primitives>> ] bi "Up to %d: %d triples, %d primitives.\n" printf ;
- pyth ( -- )
8 [1,b] [ 10^ dup count-triplets pprint-triplet-count ] each ;</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. Running time: 57.968821207 seconds
Forth
<lang Forth >
\ Two methods to create Pythagorean Triples
\ this code has been tested using Win32Forth and gforth
- pythag_fibo ( f1 f0 -- )
\ Create Pythagorean Triples from 4 element Fibonacci series
\ this is called with the first two members of a 4 element Fibonacci series
\ Price and Burkhart have two good articles about this method
\ "Pythagorean Tree: A New Species" and
\ "Heron's Formula, Descartes Circles, and Pythagorean Triangles"
\ Horadam found out how to compute Pythagorean Triples from Fibonacci series
\ compute the two other members of the Fibonacci series and put them in
\ local variables. I was unable to do this with out using locals
2DUP + 2DUP + 2OVER 2DUP + 2DUP +
LOCALS| f3 f2 f1 f0 |
wk_level @ 9 .r f0 8 .r f1 8 .r f2 8 .r f3 8 .r
\ this block calculates the sides of the Pythagorean Triangle using single precision
\ f0 f3 * 14 .r \ side a (always odd)
\ 2 f1 * f2 * 10 .r \ side b (a multiple of 4)
\ f0 f2 * f1 f3 * + 10 .r \ side c, the hyponenuse, (always odd)
\ this block calculates double precision values
f0 f3 um* 15 d.r \ side a (always odd)
2 f1 * f2 um* 15 d.r \ side b (a multiple of 4)
f0 f2 um* f1 f3 um* d+ 17 d.r cr \ side c, the hypotenuse, (always odd)
MAX_LEVEL @ wk_LEVEL @ U> IF \ TRUE if MAX_LEVEL > WK_LEVEL
wk_level @ 1+ wk_level !
\ this creates a teranary tree of Pythagorean triples
\ use a two of the members of the Fibonacci series as seeds for the
\ next level
\ It's the same tree created by Barning or Hall using matrix multiplication
f3 f1 recurse
f3 f2 recurse
f0 f2 recurse
wk_level @ 1- wk_level !
else
then
drop drop drop drop ;
\ implements the Fibonacci series -- Pythagorean triple \ the stack contents sets how many iteration levels there will be
- pf_test
\ the stack contents set up the maximum level
max_level !
0 wk_level !
cr
\ call the function with the first two elements of the base Fibonacci series
1 1 pythag_fibo ;
- gcd ( a b -- gcd )
begin ?dup while tuck mod repeat ;
\ this is the classical algorithm, known to Euclid, it is explained in many \ books on Number Theory \ this generates all primitive Pythagorean triples
\ i -- inner loop index or current loop index \ j -- outer loop index \ stack contents is the upper limit for j \ i and j can not both be odd \ the gcd( i, j ) must be 1 \ j is greater than i \ the stack contains the upper limit of the j variable
- pythag_ancn ( limit -- )
cr
1 + 2 do
i 1 and if 2 else 1 then
\ this sets the start value of the inner loop so that
\ if the outer loop index is odd only even inner loop indices happen
\ if the outer loop index is even only odd inner loop indices happen
i swap do
i j gcd 1 - 0> if else \ do this if gcd( i, j ) is 1
j 5 .r i 5 .r
\ j j * i i * - 12 .r \ a side of Pythagorean triangle (always odd)
\ i j * 2 * 9 .r \ b side of Pythagorean triangle (multiple of 4)
\ i i * j j * + 9 .r \ hypotenuse of Pythagorean triangle (always odd)
\ this block calculates double precision Pythagorean triple values
j j um* i i um* d- 15 d.r \ a side of Pythagorean triangle (always odd)
i j um* d2* 15 d.r \ b side of Pythagorean triangle (multiple of 4)
i i um* j j um* d+ 17 d.r \ hypotenuse of Pythagorean triangle (always odd)
cr then 2 +loop \ keep i being all odd or all even
loop ;
Current directory: C:\Forth ok FLOAD 'C:\Forth\ancien_fibo_pythag.F' ok
ok
ok ok
3 pf_test
0 1 1 2 3 3 4 5
1 3 1 4 5 15 8 17
2 5 1 6 7 35 12 37
3 7 1 8 9 63 16 65
3 7 6 13 19 133 156 205
3 5 6 11 17 85 132 157
2 5 4 9 13 65 72 97
3 13 4 17 21 273 136 305
3 13 9 22 31 403 396 565
3 5 9 14 23 115 252 277
2 3 4 7 11 33 56 65
3 11 4 15 19 209 120 241
3 11 7 18 25 275 252 373
3 3 7 10 17 51 140 149
1 3 2 5 7 21 20 29
2 7 2 9 11 77 36 85
3 11 2 13 15 165 52 173
3 11 9 20 29 319 360 481
3 7 9 16 25 175 288 337
2 7 5 12 17 119 120 169
3 17 5 22 27 459 220 509
3 17 12 29 41 697 696 985
3 7 12 19 31 217 456 505
2 3 5 8 13 39 80 89
3 13 5 18 23 299 180 349
3 13 8 21 29 377 336 505
3 3 8 11 19 57 176 185
1 1 2 3 5 5 12 13
2 5 2 7 9 45 28 53
3 9 2 11 13 117 44 125
3 9 7 16 23 207 224 305
3 5 7 12 19 95 168 193
2 5 3 8 11 55 48 73
3 11 3 14 17 187 84 205
3 11 8 19 27 297 304 425
3 5 8 13 21 105 208 233
2 1 3 4 7 7 24 25
3 7 3 10 13 91 60 109
3 7 4 11 15 105 88 137
3 1 4 5 9 9 40 41
ok
ok
10 pythag_ancn
2 1 3 4 5 3 2 5 12 13 4 1 15 8 17 4 3 7 24 25 5 2 21 20 29 5 4 9 40 41 6 1 35 12 37 6 5 11 60 61 7 2 45 28 53 7 4 33 56 65 7 6 13 84 85 8 1 63 16 65 8 3 55 48 73 8 5 39 80 89 8 7 15 112 113 9 2 77 36 85 9 4 65 72 97 9 8 17 144 145 10 1 99 20 101 10 3 91 60 109 10 7 51 140 149 10 9 19 180 181 ok
</lang>
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
Go
<lang go>package main
import "fmt"
var total, prim, maxPeri int64
func newTri(s0, s1, s2 int64) {
if p := s0 + s1 + s2; p <= maxPeri {
prim++
total += maxPeri / p
newTri(+1*s0-2*s1+2*s2, +2*s0-1*s1+2*s2, +2*s0-2*s1+3*s2)
newTri(+1*s0+2*s1+2*s2, +2*s0+1*s1+2*s2, +2*s0+2*s1+3*s2)
newTri(-1*s0+2*s1+2*s2, -2*s0+1*s1+2*s2, -2*s0+2*s1+3*s2)
}
}
func main() {
for maxPeri = 100; maxPeri <= 1e11; maxPeri *= 10 {
prim = 0
total = 0
newTri(3, 4, 5)
fmt.Printf("Up to %d: %d triples, %d primitives\n",
maxPeri, total, prim)
}
}</lang> Output:
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 Up to 100000000000: 161750315680 triples, 7023049293 primitives
Groovy
Parent/Child Algorithm
Solution: <lang groovy>class Triple {
BigInteger a, b, c
def getPerimeter() { this.with { a + b + c } }
boolean isValid() { this.with { a*a + b*b == c*c } }
}
def initCounts (def n = 10) {
(n..1).collect { 10g**it }.inject ([:]) { Map map, BigInteger perimeterLimit ->
map << [(perimeterLimit): [primative: 0g, total: 0g]]
}
}
def findPythagTriples, findChildTriples
findPythagTriples = {Triple t = new Triple(a:3, b:4, c:5), Map counts = initCounts() ->
def p = t.perimeter
def currentCounts = counts.findAll { pLimit, tripleCounts -> p <= pLimit }
if (! currentCounts || ! t.valid) { return }
currentCounts.each { pLimit, tripleCounts ->
tripleCounts.with { primative ++; total += pLimit.intdiv(p) }
}
findChildTriples(t, currentCounts)
counts
}
findChildTriples = { Triple t, Map counts ->
t.with {
[
[ a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c],
[ a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c],
[-a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c]
]*.sort().each { aa, bb, cc ->
findPythagTriples(new Triple(a:aa, b:bb, c:cc), counts)
}
}
}</lang>
Test: <lang groovy>printf (' LIMIT PRIMATIVE ALL\n') findPythagTriples().sort().each { perimeterLimit, result ->
def exponent = perimeterLimit.toString().size() - 1
printf ('a+b+c <= 10E%2d %9d %12d\n', exponent, result.primative, result.total)
}</lang>
Output:
LIMIT PRIMATIVE ALL a+b+c <= 10E 1 0 0 a+b+c <= 10E 2 7 17 a+b+c <= 10E 3 70 325 a+b+c <= 10E 4 703 4858 a+b+c <= 10E 5 7026 64741 a+b+c <= 10E 6 70229 808950 a+b+c <= 10E 7 702309 9706567 a+b+c <= 10E 8 7023027 113236940 a+b+c <= 10E 9 70230484 1294080089 a+b+c <= 10E10 702304875 14557915466
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
Haskell
<lang haskell>pytr n = filter (\(_, a, b, c) -> a+b+c <= n) [(prim a b c, a, b, c) | a <- [1..n], b <- [1..n], c <- [1..n], a < b && b < c, a^2 + b^2 == c^2]
where prim a b _ = gcd a b == 1
main = putStrLn $ "Up to 100 there are " ++ (show $ length xs) ++ " triples, of which " ++ (show $ length $ filter (\(x,_,_,_) -> x == True) xs) ++ " are primitive."
where xs = pytr 100
</lang>
Up to 100 there are 17 triples, of which 7 are primitive.
Recursive primitive generation: <lang haskell>triangles :: Int -> Int triangles max_peri | max_peri < 12 = []
| otherwise = concat tiers where tiers = takeWhile (not.null) $ iterate tier 3,4,5 tier = concatMap (filter ((<=max_peri).sum).tmul) tmul t = map (map (sum . zipWith (*) t)) [[[ 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]]]
triangle_count max_p = (length t, sum $ map ((max_p `div`).sum) t)
where t = triangles max_p
main = mapM_ putStrLn $
map (\n -> show n ++ " " ++ show (triangle_count n)) $ map (10^) [1..]</lang>
- Output:
10 (0,0) 100 (7,17) 1000 (70,325) 10000 (703,4858) 100000 (7026,64741) 1000000 (70229,808950) 10000000 (702309,9706567) ...
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),
a2.negate().add(b).add(c2),
a2.negate().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.
jq
The jq program presented here is based on Euclid's formula, and uses the same algorithm and notation as in the AutoHotKey section.
The implementation illustrates how an inner function with arity 0 can attain a high level of efficiency with both jq 1.4 and later. A simpler implementation is possible with versions of jq greater than 1.4. <lang jq>def gcd(a; b):
def _gcd: if .[1] == 0 then .[0] else [.[1], .[0] % .[1]] | _gcd end; [a,b] | _gcd ;
- Return: [total, primitives] for pythagorean triangles having
- perimeter no larger than peri.
- The following uses Euclid's formula with the convention: m > n.
def count(peri):
# The inner function can be used to count for a given value of m:
def _count:
# state [n,m,p, [total, primitives]]
.[0] as $n | .[1] as $m | .[2] as $p
| if $n < $m and $p <= peri then
if (gcd($m;$n) == 1)
then .[3] | [ (.[0] + ((peri/$p)|floor) ), (.[1] + 1)]
else .[3]
end
| [$n+2, $m, $p+4*$m, .] | _count
else .
end;
# m^2 < m*(m+1) <= m*(m+n) = perimeter/2
reduce range(2; (peri/2) | sqrt + 1) as $m
( [1, 2, 12, [0,0]];
(($m % 2) + 1) as $n
| (2 * $m * ($m + $n) ) as $p # a+b+c for this (m,n)
| [$n, $m, $p, .[3]] | _count
) | .[3] ;
- Example:
<lang jq>def pow(i): . as $in | reduce range(0; i) as $j (1; . * $in);
range(1; 9) | . as $i | 10|pow($i) as $i | "\($i): \(count($i) )" </lang>
- Output:
<lang sh>$ jq -M -c -r -n -f Pythagorean_triples.jq 10: [0,0] 100: [17,7] 1000: [325,70] 10000: [4858,703] 100000: [64741,7026] 1000000: [808950,70229] 10000000: [9706567,702309] 100000000: [113236940,7023027] </lang>
Julia
This solution uses the the Euclidian concept of m and n as generators of Pythagorean triplets. When m and n are coprime and have opposite parity, the generated triplets are primitive. It works reasonably well up to a limit of 10^10. <lang Julia> function primitiven{T<:Integer}(m::T)
1 < m || return T[]
m != 2 || return T[1]
!isprime(m) || return T[2:2:m-1]
rp = trues(m-1)
if isodd(m)
rp[1:2:m-1] = false
end
for p in keys(factor(m))
rp[p:p:m-1] = false
end
T[1:m-1][rp]
end
function pythagoreantripcount{T<:Integer}(plim::T)
primcnt = 0
fullcnt = 0
11 < plim || return (primcnt, fullcnt)
for m in 2:plim
p = 2m^2
p+2m <= plim || break
for n in primitiven(m)
q = p + 2m*n
q <= plim || break
primcnt += 1
fullcnt += div(plim, q)
end
end
return (primcnt, fullcnt)
end
println("Counting Pythagorian Triplets within perimeter limits:") println(" Limit All Primitive") for om in 1:10
(pcnt, fcnt) = pythagoreantripcount(10^om) println(@sprintf " 10^%02d %11d %9d" om fcnt pcnt)
end </lang>
- Output:
Counting Pythagorian Triplets within perimeter limits:
Limit All Primitive
10^01 0 0
10^02 17 7
10^03 325 70
10^04 4858 703
10^05 64741 7026
10^06 808950 70229
10^07 9706567 702309
10^08 113236940 7023027
10^09 1294080089 70230484
10^10 14557915466 702304875
Lasso
<lang lasso>// Brute Force: Too slow for large numbers define num_pythagorean_triples(max_perimeter::integer) => {
local(max_b) = (#max_perimeter / 3)*2
return (
with a in 1 to (#max_b - 1)
sum integer(
with b in (#a + 1) to #max_b
let c = math_sqrt(#a*#a + #b*#b)
where #c == integer(#c)
where #c > #b
where (#a+#b+#c) <= #max_perimeter
sum 1
)
)
} stdout(`Number of Pythagorean Triples in a Perimeter of 100: `) stdoutnl(num_pythagorean_triples(100)) </lang> Output:
Number of Pythagorean Triples in a Perimeter of 100: 17
Liberty BASIC
<lang lb> print time$()
for power =1 to 6
perimeterLimit =10^power upperBound =int( 1 +perimeterLimit^0.5) primitives =0 triples =0 extras =0 ' will count the in-range multiples of any primitive
for m =2 to upperBound
for n =1 +( m mod 2 =1) to m -1 step 2
term1 =2 *m *n
term2 =m *m -n *n
term3 =m *m +n *n
perimeter =term1 +term2 +term3
if perimeter <=perimeterLimit then triples =triples +1
a =term1
b =term2
do
r = a mod b
a =b
b =r
loop until r <=0
if ( a =1) and ( perimeter <=perimeterLimit) then 'we've found a primitive triple if a =1, since hcf =1. And it is inside perimeter range. Save it in an array
primitives =primitives +1
if term1 >term2 then temp =term1: term1 =term2: term2 =temp 'swap so in increasing order of side length
nEx =int( perimeterLimit /perimeter) 'We have the primitive & removed any multiples. Now calculate ALL the multiples in range.
extras =extras +nEx
end if
scan
next n
next m
print " Number of primitives having perimeter below "; 10^power, " was "; primitives, " & "; extras, " non-primitive triples." print time$()
next power
print "End" end </lang>
17:59:34 Number of primitives having perimeter below 10 was 0 & 0 non-primitive triples. 17:59:34 Number of primitives having perimeter below 100 was 7 & 17 non-primitive triples. 17:59:34 Number of primitives having perimeter below 1000 was 70 & 325 non-primitive triples. 17:59:34 Number of primitives having perimeter below 10000 was 703 & 4858 non-primitive triples. 17:59:35 Number of primitives having perimeter below 100000 was 7026 & 64741 non-primitive triples. 17:59:42 Number of primitives having perimeter below 1000000 was 70229 & 808950 non-primitive triples. 18:01:30 End
Mathematica
Short code but not a very scalable approach... <lang Mathematica>pythag[n_] := Block[{soln = Solve[{a^2 + b^2 == c^2, a + b + c <= n, 0 < a < b < c}, {a, b, c}, Integers]},
{Length[soln], Count[GCD[a, b] == GCD[b, c] == GCD[c, a] == 1 /. soln, True]}
]</lang>
- Output:
pythag[10]
{0,0}
pythag[100]
{17, 7}
pythag[1000]
{325, 70}
MATLAB / Octave
<lang Matlab>N= 100;
a = 1:N; b = a(ones(N,1),:).^2; b = b+b'; b = sqrt(b); [y,x]=find(b==fix(b)); % test % here some alternative tests % b = b.^(1/k); [y,x]=find(b==fix(b)); % test 2 % [y,x]=find(b==(fix(b.^(1/k)).^k)); % test 3 % b=b.^(1/k); [y,x]=find(abs(b - round(b)) <= 4*eps*b); z = sqrt(x.^2+y.^2); ix = (z+x+y<100) & (x < y) & (y < z); p = find(gcd(x(ix),y(ix))==1); % find primitive triples
printf('There are %i Pythagorean Triples and %i primitive triples with a perimeter smaller than %i.\n',...
sum(ix), length(p), N); </lang>
Output:
There are 17 Pythagorean Triples and 7 primitive triples with a perimeter smaller than 100.
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
Nim
<lang nim>const u = [[ 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]]
var
total, prim = 0 maxPeri = 10
proc newTri(ins) =
var p = ins[0] + ins[1] + ins[2] if p > maxPeri: return inc(prim) total += maxPeri div p
for i in 0..2:
newTri([u[i][0] * ins[0] + u[i][1] * ins[1] + u[i][2] * ins[2],
u[i][3] * ins[0] + u[i][4] * ins[1] + u[i][5] * ins[2],
u[i][6] * ins[0] + u[i][7] * ins[1] + u[i][8] * ins[2]])
while maxPeri <= 100_000_000:
total = 0 prim = 0 newTri([3, 4, 5]) echo "Up to ", maxPeri, ": ", total, " triples, ", prim, " primitives" maxPeri *= 10</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
OCaml
<lang OCaml>let isqrt n =
let rec iter t =
let d = n - t*t in
if (0 <= d) && (d < t+t+1) (* t*t <= n < (t+1)*(t+1) *)
then t else iter ((t+(n/t))/2)
in iter 1
let rec gcd a b =
let t = a mod b in if t = 0 then b else gcd b t
let coprime a b = gcd a b = 1
let num_to ms =
let ctr = ref 0 in
let prim_ctr = ref 0 in
let max_m = isqrt (ms/2) in
for m = 2 to max_m do
for j = 0 to (m/2) - 1 do
let n = m-(2*j+1) in
if coprime m n then
let s = 2*m*(m+n) in
if s <= ms then
(ctr := !ctr + (ms/s); incr prim_ctr)
done
done;
(!ctr, !prim_ctr)
let show i =
let s, p = num_to i in Printf.printf "For perimeters up to %d there are %d total and %d primitive\n%!" i s p;;
List.iter show [ 100; 1000; 10000; 100000; 1000000; 10000000; 100000000 ]</lang> Output:
For perimeters up to 100 there are 17 total and 7 primitive For perimeters up to 1000 there are 325 total and 70 primitive For perimeters up to 10000 there are 4858 total and 703 primitive For perimeters up to 100000 there are 64741 total and 7026 primitive For perimeters up to 1000000 there are 808950 total and 70229 primitive For perimeters up to 10000000 there are 9706567 total and 702309 primitive For perimeters up to 100000000 there are 113236940 total and 7023027 primitive
PARI/GP
This version is reasonably efficient and can handle inputs like a million quickly. <lang parigp>do(lim)={
my(prim,total,P);
lim\=1;
for(m=2,sqrtint(lim\2),
forstep(n=1+m%2,min(sqrtint(lim-m^2),m-1),2,
P=2*m*(m+n);
if(gcd(m,n)==1 && P<=lim,
prim++;
total+=lim\P
)
)
);
[prim,total]
}; do(100)</lang>
Pascal
<lang pascal>Program PythagoreanTriples (output);
var
total, prim, maxPeri: int64;
procedure newTri(s0, s1, s2: int64);
var
p: int64;
begin
p := s0 + s1 + s2;
if p <= maxPeri then
begin
inc(prim);
total := total + maxPeri div p;
newTri( s0 + 2*(-s1+s2), 2*( s0+s2) - s1, 2*( s0-s1+s2) + s2);
newTri( s0 + 2*( s1+s2), 2*( s0+s2) + s1, 2*( s0+s1+s2) + s2);
newTri(-s0 + 2*( s1+s2), 2*(-s0+s2) + s1, 2*(-s0+s1+s2) + s2);
end;
end;
begin
maxPeri := 100;
while maxPeri <= 1e10 do
begin
prim := 0;
total := 0;
newTri(3, 4, 5);
writeln('Up to ', maxPeri, ': ', total, ' triples, ', prim, ' primitives.');
maxPeri := maxPeri * 10;
end;
end.</lang> Output (on Core2Duo 2GHz laptop):
time ./PythagoreanTriples 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. 109.694u 0.094s 1:50.43 99.4% 0+0k 0+0io 0pf+0w
Perl
<lang perl>sub gcd {
my ($n, $m) = @_;
while($n){
my $t = $n;
$n = $m % $n;
$m = $t;
}
return $m;
}
sub tripel {
my $pmax = shift;
my $prim = 0;
my $count = 0;
my $nmax = sqrt($pmax)/2;
for( my $n=1; $n<=$nmax; $n++ ) {
for( my $m=$n+1; (my $p = 2*$m*($m+$n)) <= $pmax; $m+=2 ) {
next unless 1==gcd($m,$n);
$prim++;
$count += int $pmax/$p;
}
}
printf "Max. perimeter: %d, Total: %d, Primitive: %d\n", $pmax, $count, $prim;
}
tripel 10**$_ for 1..8; </lang>
- Output:
Max. perimeter: 10, Total: 0, Primitive: 0 Max. perimeter: 100, Total: 17, Primitive: 7 Max. perimeter: 1000, Total: 325, Primitive: 70 Max. perimeter: 10000, Total: 4858, Primitive: 703 Max. perimeter: 100000, Total: 64741, Primitive: 7026 Max. perimeter: 1000000, Total: 808950, Primitive: 70229 Max. perimeter: 10000000, Total: 9706567, Primitive: 702309 Max. perimeter: 100000000, Total: 113236940, Primitive: 7023027
Perl 6
Here is a straight-forward, naive brute force implementation: <lang perl6>constant limit = 100;
for [X] [^limit] xx 3 -> (\a, \b, \c) {
say [a, b, c] if a < b < c and a**2 + b**2 == c**2
}</lang>
Here is a slightly less naive brute force implementation that is not really practical for large perimeter limits. It's pretty zippy up to about 10000 though.
<lang perl6>my %triples; my $limit = 10000;
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>
- Output:
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 ... 196 4800 4804 310 4800 4810 392 4794 4810 171 4872 4875 99 4900 4901 - primitive 140 4899 4901 - primitive There are 4858 Pythagorean Triples with a perimeter <= 10000, of which 703 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.
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.
PL/I
Version 1 <lang PL/I>*process source attributes xref or(!);
/*********************************************************************
* REXX pgm counts number of Pythagorean triples
* that exist given a max perimeter of N,
* and also counts how many of them are primatives.
* 05.05.2013 Walter Pachl translated from REXX version 2
*********************************************************************/
pyt: Proc Options(main);
Dcl sysprint Print;
Dcl (addr,mod,right) Builtin;
Dcl memn Bin Fixed(31) Init(0);
Dcl mabca(300) Char(12);
Dcl 1 mabc,
2 ma Dec fixed(7),
2 mb Dec fixed(7),
2 mc Dec fixed(7);
Dcl mabce Char(12) Based(addr(mabc));
Dcl 1 abc,
2 a Dec fixed(7),
2 b Dec fixed(7),
2 c Dec fixed(7);
Dcl abce Char(12) Based(addr(abc));
Dcl (prims,trips,m,n,aa,aabb,cc,aeven,ab) Dec Fixed(7);
mabca=;
trips=0;
prims=0;
n=100;
la:
Do a=3 To n/3;
aa=a*a; /* limit side to 1/3 of perimeter.*/
aeven=mod(a,2)=0;
lb:Do b=a+1 By 1+aeven; /* triangle can't be isosceles. */
ab=a+b; /* compute partial perimeter. */
If ab>=n Then
Iterate la; /* a+b>perimeter? Try different A*/
aabb=aa+b*b; /* compute sum of a² + b² (cheat)*/
Do c=b+1 By 1;
cc=c*c; /* 3rd side: also compute c² */
If aeven Then
If mod(c,2)=0 Then
Iterate;
If ab+c>n Then
Iterate la; /* a+b+c > perimeter? Try diff A.*/
If cc>aabb Then
Iterate lb; /* c² > a²+b² ? Try different B.*/
If cc^=aabb Then
Iterate; /* c² ¬= a²+b² ? Try different C.*/
If mema(abce) Then
Iterate;
trips=trips+1; /* eureka. */
prims=prims+1; /* count this primitive triple. */
Put Edit(a,b,c,' ',right(a**2+b**2,5),right(c**2,5),a+b+c)
(Skip,f(4),2(f(5)),a,2(f(6)),f(9));
Do m=2 By 1;
ma=a*m;
mb=b*m;
mc=c*m; /* gen non-primitives. */
If ma+mb+mc>n Then
Leave;
/* is this multiple a triple ? */
trips=trips+1; /* yuppers, then we found another.*/
If mod(m,2)=1 Then /* store as even multiple. */
call mems(mabce);
Put Edit(ma,mb,mc,' * ',
right(ma**2+mb**2,5),right(mc**2,5),ma+mb+mc)
(Skip,f(4),2(f(5)),a,2(f(6)),f(9));
End; /* m */
End; /* c */
End; /* b */
End; /* a */
Put Edit('max perimeter = ',n, /* show a single line of output. */
'Pythagorean triples =',trips,
'primitives =',prims)
(Skip,a,f(5),2(x(9),a,f(4)));
mems: Proc(e); Dcl e Char(12); memn+=1; mabca(memn)=e; End;
mema: Proc(e) Returns(bit(1));
Dcl e Char(12);
Do memi=1 To memn;
If mabca(memi)=e Then Return('1'b);
End;
Return('0'b);
End;
End;</lang>
Version 2 <lang PL/I> pythagorean: procedure options (main, reorder); /* 23 January 2014 */
declare (a, b, c) fixed (3); declare (asquared, bsquared) fixed; declare (triples, primitives) fixed binary(31) initial (0);
do a = 1 to 100;
asquared = a*a;
do b = a+1 to 100;
bsquared = b*b;
do c = b+1 to 100;
if a+b+c <= 100 then
if asquared + bsquared = c*c then
do;
triples = triples + 1;
if GCD(a,b) = 1 then primitives = primitives + 1;
end;
end;
end;
end;
put skip data (triples, primitives);
GCD: procedure (a, b) returns (fixed binary (31)) recursive;
declare (a, b) fixed binary (31);
if b = 0 then return (a);
return (GCD (b, mod(a, b)) );
end GCD;
end pythagorean;</lang> Output:
TRIPLES= 17 PRIMITIVES= 7;
Output for P=1000:
TRIPLES= 325 PRIMITIVES= 70;
PowerShell
<lang PowerShell> function triples($p) {
if($p -gt 4) {
# ai + bi + ci = pi <= p
# ai < bi < ci --> 3ai < pi <= p and ai + 2bi < pi <= p
$pa = [Math]::Floor($p/3)
1..$pa | foreach {
$ai = $_
$pb = [Math]::Floor(($p-$ai)/2)
($ai+1)..$pb | foreach {
$bi = $_
$pc = $p-$ai-$bi
($bi+1)..$pc | where {
$ci = $_
$pi = $ai + $bi + $ci
$ci*$ci -eq $ai*$ai + $bi*$bi
} |
foreach {
[pscustomobject]@{
a = "$ai"
b = "$bi"
c = "$ci"
p = "$pi"
}
}
}
}
}
else {
Write-Error "$p is not greater than 4"
}
} function gcd ($a, $b) {
function pgcd ($n, $m) {
if($n -le $m) {
if($n -eq 0) {$m}
else{pgcd $n ($m%$n)}
}
else {pgcd $m $n}
}
$n = [Math]::Abs($a)
$m = [Math]::Abs($b)
(pgcd $n $m)
} $triples = (triples 100)
$coprime = $triples | where {((gcd $_.a $_.b) -eq 1) -and ((gcd $_.a $_.c) -eq 1) -and ((gcd $_.b $_.c) -eq 1)}
"There are $(($triples).Count) Pythagorean triples with perimeter no larger than 100
and $(($coprime).Count) of them are coprime."
</lang> Output:
There are 17 Pythagorean triples with perimeter no larger than 100 and 7 of them are coprime.
PureBasic
<lang purebasic>
Procedure.i ConsoleWrite(t.s) ; compile using /CONSOLE option
OpenConsole()
PrintN (t.s)
CloseConsole()
ProcedureReturn 1
EndProcedure
Procedure.i StdOut(t.s) ; compile using /CONSOLE option
OpenConsole()
Print(t.s)
CloseConsole()
ProcedureReturn 1
EndProcedure
Procedure.i gcDiv(n,m) ; greatest common divisor if n=0:ProcedureReturn m:endif while m <> 0
if n > m n - m else m - n endif
wend ProcedureReturn n EndProcedure
st=ElapsedMilliseconds()
nmax =10000 power =8
dim primitiveA(power) dim alltripleA(power) dim pmaxA(power)
x=1 for i=1 to power
x*10:pmaxA(i)=x/2
next
for n=1 to nmax
for m=(n+1) to (nmax+1) step 2 ; assure m-n is odd
d=gcDiv(n,m)
p=m*m+m*n
for i=1 to power
if p<=pmaxA(i)
if d =1
primitiveA(i)+1 ; right here we have the primitive perimeter "seed" 'p'
k=1:q=p*k ; set k to one to include p : use q as the 'p*k'
while q<=pmaxA(i)
alltripleA(i)+1 ; accumulate multiples of this perimeter while q <= pmaxA(i)
k+1:q=p*k
wend
endif
endif
next
next
next
for i=1 to power
t.s="Up to "+str(pmaxA(i)*2)+": " t.s+str(alltripleA(i))+" triples, " t.s+str(primitiveA(i))+" primitives." ConsoleWrite(t.s)
next ConsoleWrite("") et=ElapsedMilliseconds()-st:ConsoleWrite("Elapsed time = "+str(et)+" milliseconds") </lang>
- Output
>cmd /c "C:\_sys\temp\PythagoreanTriples.exe" 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. Elapsed time = 5163 milliseconds >Exit code: 0
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>
Racket
<lang racket>#lang racket
- | Euclid's enumeration formula and counting is fast enough for extra credit.
For maximum perimeter P₀, the primitive triples are enumerated by n,m with: 1 ≤ n < m perimeter P(n, m) ≤ P₀ where P(n, m) = (m² - n²) + 2mn + (m² + n²) = 2m(m+n) m and n of different parity and coprime.
Since n < m, a simple close non-tight bound on n is P(n, n) < P₀. For each of these the exact set of m's can be enumerated.
Each primitive triple with perimeter p represents one triple for each kp ≤ P₀, of which there are floor(P₀/p) k's. |#
(define (P n m) (* 2 m (+ m n))) (define (number-of-triples P₀)
(for/fold ([primitive 0] [all 0])
([n (in-naturals 1)]
#:break (>= (P n n) P₀))
(for*/fold ([primitive′ primitive] [all′ all])
([m (in-naturals (+ n 1))]
#:break (> (P n m) P₀)
#:when (and (odd? (- m n)) (coprime? m n)))
(values (+ primitive′ 1)
(+ all′ (quotient P₀ (P n m)))))))
(define (print-results P₀)
(define-values (primitive all) (number-of-triples P₀))
(printf "~a ~a:\n ~a, ~a.\n"
"Number of Pythagorean triples and primitive triples with perimeter ≤"
P₀
all primitive))
(print-results 100) (time (print-results (* 100 1000 1000)))
- |
Number of Pythagorean triples and primitive triples with perimeter ≤ 100:
17, 7.
Number of Pythagorean triples and primitive triples with perimeter ≤ 100000000:
113236940, 7023027.
cpu time: 11976 real time: 12215 gc time: 2381
|#</lang>
REXX
using GCD for determinacy
<lang rexx>/*REXX program counts number of Pythagorean triples that exist given a max */ /*────────── perimeter of N, and also counts how many of them are primitives.*/ trips=0; prims=0 /*set the number of triples, primitives*/ parse arg N .; if N== then n=100 /*N specified? No, then use default. */
do a=3 to N%3; aa=a*a /*limit side to 1/3 of the perimeter.*/
do b=a+1 /*the triangle can't be isosceles. */
ab=a+b /*compute a partial perimeter (2 sides)*/
if ab>=N then iterate a /*is a+b ≥ perimeter? Try different A*/
aabb=aa+b*b /*compute the sum of a²+b² (shortcut)*/
do c=b+1 /*compute the value of the third side. */
if ab+c>N then iterate a /*is a+b+c > perimeter? Try diff. A.*/
cc=c*c /*compute the value of C². */
if cc > aabb then iterate b /*is c² > a²+b² ? Try a different B.*/
if cc\==aabb then iterate /*is c² ¬= a²+b² ? Try a different C.*/
trips=trips+1 /*eureka. We found a Pythagorean triple*/
prims=prims+(gcd(a,b)==1) /*is this triple a primitive triple? */
end /*c*/
end /*b*/
end /*a*/
_=left(,7) /*for padding the output with 7 blanks.*/ say 'max perimeter =' N _ "Pythagorean triples =" trips _ 'primitives =' prims exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────────────*/ gcd: procedure; parse arg x,y; do until y==0; parse value x//y y with y x; end; return x</lang> output when using the default input of: 100
max perimeter = 100 Pythagorean triples = 17 primitives = 7
output when using the input of: 1000
max perimeter = 1000 Pythagorean triples = 325 primitives = 70
using single evenness for determinacy
This REXX version takes advantage that primitive Pythagorean triples must have one and only one even number.
Non-primitive Pythagorean triples are generated after a primitive triple is found.
<lang rexx>/*REXX program counts number of Pythagorean triples that exist given a max */
/*REXX program counts number of Pythagorean triples that exist given a max */
/*────────── perimeter of N, and also counts how many of them are primitives.*/
@.=0; trips=0; prims=0 /*define some REXX variables to zero. */
parse arg N .; if N== then n=100 /*N specified? No, then use default. */
do a=3 to N%3; aa=a*a /*limit side to 1/3 of the perimeter.*/ aEven= a//2==0 /*set variable to 1 if A is even. */
do b=a+1 by 1+aEven /*the triangle can't be isosceles. */
ab=a+b /*compute a partial perimeter (2 sides)*/
if ab>=N then iterate a /*is a+b ≥ perimeter? Try different A*/
aabb=aa+b*b /*compute the sum of a²+b² (shortcut)*/
do c=b+1 /*compute the value of the third side. */
if aEven then if c//2==0 then iterate
if ab+c>n then iterate a /*a+b+c > perimeter? Try different A.*/
cc=c*c /*compute the value of C². */
if cc > aabb then iterate b /*is c² > a²+b² ? Try a different B.*/
if cc\==aabb then iterate /*is c² ¬= a²+b² ? Try a different C.*/
if @.a.b.c then iterate /*Is this a duplicate? Then try again.*/
trips=trips+1 /*Eureka! We found a Pythagorean triple*/
prims=prims+1 /*count this also as a primitive triple*/
do m=2; am=a*m; bm=b*m; cm=c*m /*generate non-primitives.*/
if am+bm+cm>N then leave /*is this multiple Pythagorean triple? */
trips=trips+1 /*Eureka! We found a Pythagorean triple*/
@.am.bm.cm=1 /*mark Pythagorean triangle as a triple*/
end /*m*/
end /*c*/
end /*b*/
end /*a*/
_=left(,7) /*for padding the output with 7 blanks.*/ say 'max perimeter =' N _ "Pythagorean triples =" trips _ 'primitives =' prims
/*stick a fork in it, we're all done. */</lang>
output is identical to the 1st REXX version.
Ruby
<lang ruby>class PythagoranTriplesCounter
def initialize(limit) @limit = limit @total = 0 @primitives = 0 generate_triples(3, 4, 5) end attr_reader :total, :primitives private def generate_triples(a, b, c) perim = a + b + c return if perim > @limit
@primitives += 1 @total += @limit / perim
generate_triples( a-2*b+2*c, 2*a-b+2*c, 2*a-2*b+3*c) generate_triples( a+2*b+2*c, 2*a+b+2*c, 2*a+2*b+3*c) generate_triples(-a+2*b+2*c,-2*a+b+2*c,-2*a+2*b+3*c) end
end
perim = 10 while perim <= 100_000_000
c = PythagoranTriplesCounter.new perim p [perim, c.total, c.primitives] perim *= 10
end</lang>
output
[10, 0, 0] [100, 17, 7] [1000, 325, 70] [10000, 4858, 703] [100000, 64741, 7026] [1000000, 808950, 70229] [10000000, 9706567, 702309] [100000000, 113236940, 7023027]
Scheme
<lang Scheme>(use srfi-42)
(define (py perim)
(define prim 0)
(values
(sum-ec
(: c perim) (: b c) (: a b)
(if (and (<= (+ a b c) perim)
(= (square c) (+ (square b) (square a)))))
(begin (when (= 1 (gcd a b)) (inc! prim)))
1)
prim))</lang>
Testing:
gosh> (py 100) 17 7
Scratch
Scratch is a visual programming language. Click the link, then "see inside" to see the code.
https://scratch.mit.edu/projects/79066598/
Output: 17 Pythagorean triples with a perimeter less than 100, 7 of which are primitive.
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.
Sidef
<lang ruby>func triples(limit) {
var primitive = 0; var civilized = 0;
func oyako(a, b, c) {
(var perim = a+b+c) > limit || (
primitive++;
civilized += int(limit / 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})";
}
range(1, Math.inf).each { |n|
say triples(10**n);
};</lang>
- Output:
10 => (0 0) 100 => (7 17) 1000 => (70 325) 10000 => (703 4858) 100000 => (7026 64741) ^C
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
VBScript
<lang vb> For i=1 To 8 WScript.StdOut.WriteLine triples(10^i) Next
Function triples(pmax) prim=0 : count=0 : nmax=Sqr(pmax)/2 : n=1 Do While n <= nmax m=n+1 : p=2*m*(m+n) Do While p <= pmax If gcd(m,n)=1 Then prim=prim+1 count=count+Int(pmax/p) End If m=m+2 p=2*m*(m+n) Loop n=n+1 Loop triples = "Max Perimeter: " & pmax &_ ", Total: " & count &_ ", Primitive: " & prim End Function
Function gcd(a,b) c = a : d = b Do If c Mod d > 0 Then e = c Mod d c = d d = e Else gcd = d Exit Do End If Loop End Function </lang>
zkl
<lang zkl>fcn tri(lim,a=3,b=4,c=5){
p:=a + b + c;
if(p>lim) return(0,0);
T(1,lim/p).zipWith('+,
tri(lim, a - 2*b + 2*c, 2*a - b + 2*c, 2*a - 2*b + 3*c),
tri(lim, a + 2*b + 2*c, 2*a + b + 2*c, 2*a + 2*b + 3*c),
tri(lim, -a + 2*b + 2*c, -2*a + b + 2*c, -2*a + 2*b + 3*c)
);
}</lang> <lang zkl>n:=10; do(10){ println("%,d: %s".fmt(n,tri(n).reverse())); n*=10; }</lang>
- Output:
10: L(0,0) 100: L(17,7) 1,000: L(325,70) 10,000: L(4858,703) 100,000: L(64741,7026) 1,000,000: L(808950,70229) 10,000,000: L(9706567,702309) VM#1 caught this unhandled exception: AssertionError : That is one big stack, infinite recursion? Stack trace for VM#1 (): jj.tri addr:112 args(4) reg(1) R <repeats 3578 times> jj.__constructor addr:23 args(0) reg(2) R ...
Max stack size is arbitrary but not adjustable.
- Programming Tasks
- Solutions by Programming Task
- Ada
- AutoHotkey
- BBC BASIC
- Bracmat
- C
- C sharp
- Clojure
- CoffeeScript
- Common Lisp
- D
- Eiffel
- Elixir
- Erlang
- ERRE
- Euphoria
- F Sharp
- Factor
- Forth
- Fortran
- Go
- Groovy
- Icon
- Unicon
- Icon Programming Library
- Haskell
- J
- Java
- Arbitrary precision
- Jq
- Julia
- Lasso
- Liberty BASIC
- Mathematica
- MATLAB
- Octave
- Modula-3
- Nim
- OCaml
- PARI/GP
- Pascal
- Perl
- Perl 6
- PicoLisp
- PL/I
- PowerShell
- PureBasic
- Python
- Racket
- REXX
- Ruby
- Scheme
- Scratch
- Seed7
- Sidef
- Tcl
- VBScript
- Zkl
- Geometry
- Pages that use a deprecated format of the math tags