# Primitive Pythagorean Triples

(related)

A Pythagorean Triple is a list (a, b, c) that satisfies the equation a2 + b2 = c2.

A Primitive Pythagorean Triple (PPT) is one where a, b, and c are all coprime (i.e., the only common divisor between the three elements is 1). For example, the (3, 4, 5) right triangle is a famous Primitive Pythagorean Triple.

## The Challenge

• Given input n, output the nth PPT. Or,
• Given input n, output the first n PPTs.

There are multiple ways to order these PPTs to form a well-ordered list, to determine which is the nth. You can choose any ordering you want, so long as you can prove (informally is fine) that your algorithm can generate every possible unique PPT. For example, your code should not output both (3,4,5) and (4,3,5) since those are duplicates of the same triple -- one or the other, please.

Similarly, whether your code is zero- or one-indexed is fine, so long as you state which you're using.

## Examples

For the examples below, I'm using one-indexing, outputting the nth PPT, and ordering by smallest c, then smallest a, then smallest b.

n | output
1 | (3, 4, 5)
2 | (5, 12, 13)
5 | (20, 21, 29)
12| (48, 55, 73)


## Rules

• The input and output can be given in any convenient format.
• In your submission, please state how your entries are ordered, and whether your entries are 0-indexed or 1-indexed.
• Your chosen ordering cannot create duplicates.
• Either a full program or a function are acceptable. If a function, you can return the output rather than printing it.
• If possible, please include a link to an online testing environment so other people can try out your code!
• Standard loopholes are forbidden.
• This is so all usual golfing rules apply, and the shortest code (in bytes) wins.
• Related – miles Jul 25 '17 at 13:29
• – miles Jul 25 '17 at 13:38
• What is the highest input we have to support? Can we assume that it fits the capabilities of our language of choice? – Mr. Xcoder Jul 25 '17 at 13:40
• @Mr.Xcoder Yes; that's a standard safe assumption, unless you're using that to exploit a loophole (e.g., the language only supports 1-bit numbers) to make the problem trivial. – AdmBorkBork Jul 25 '17 at 13:43
• I found the answer to my question: a and b must be coprime and this is enough proofwiki.org/wiki/… – edc65 Jul 26 '17 at 9:25

# Jelly, 27 25 bytes

2 bytes thanks to Jonathan Allan.

²IH;Pµ;ÆḊ
+2ḶḤ‘Œcg/ÐṂÇ€Ṣḣ


Try it online!

Outputs the first n 1-indexed triples [b, a, c], sorted by increasing b and then a.

Uses the algorithm from Wikipedia: This generates all primitive triples for all unique coprime pairs of odd integers m > n > 0.

### Explanation

+2ḶḤ‘Œcg/ÐṂÇ€Ṣḣ    Main link. Argument: n
+2                   Add 2 to n, to get enough results.
Ḷ                  Get integers [0, 1, ..., n+1].
Ḥ                 Double to get [0, 2, ..., 2n+2].
‘                Increment to get [1, 3, ..., 2n+3].
Œc              Get all ordered pairs [[1, 3], [1, 5], ..., [2n+1, 2n+3]].
g/            GCD of each pair.
ÐṂ          Grab the pairs with minimal GCD (which is 1).
Ç€        Call the helper link on each pair to get the triples.
Ṣ       Sort the triples first by a, then by b, then by c.
ḣ      Get the last n.

²IH;Pµ;ÆḊ    Helper link. Argument: pair [m, n]
²              Square to get [m², n²].
I             Increments: get [m²-n²].
H            Halve: get [(m²-n²)/2], i.e. [b].
P          Product: get mn, i.e. a.
;           Append to get [b, a].
µ         Begin a new monadic chain with argument [b, a].
ÆḊ      Get the length of the vector, i.e. c.
;        Append to get [b, a, c].

• That's a really nice explanation. Thanks! – AdmBorkBork Jul 25 '17 at 15:48
• g/Ị$Ðf -> g/ÐṂ to save two bytes (since the minimal gcd is 1 and there will always be at least one such entry). – Jonathan Allan Jul 25 '17 at 17:18 • Another byte may also be saved (although making it less efficient) by replacing +2ḶḤ‘Œc with ²Rm2Œc - scrap that it wont work for an input of 1 :( – Jonathan Allan Jul 25 '17 at 17:59 • @JonathanAllan Thanks for the minimum one. I did try a lot of 2-byte ranges, but unfortunately none were big enough. (²ḶḤ‘Œc was one of the first I thought of.) – PurkkaKoodari Jul 25 '17 at 18:02 # MATL, 36 bytes @:Ut&+wmR&fZd1Mhw/vXutnGE<]GY)t&2|h  Input is 1-based. The output order guarantees that each triple appears exactly once. The order is explained in the following. The explanation requires delving a little into how the program works. The code keeps increasing a counter k in a loop, starting at 1. For each k it generates all pairs with a = 1,...,k, b = 1,...,k, a < b, and picks those that give a Pythagorean triple with c <= k. The pairs are obtained in order of increasing b, then a. Each pair is then divided by its gcd. The resulting (possibly duplicated) pairs are arranged as a two-column matrix. This matrix is vertically concatenated with a similar matrix containing the accumulated results obtained for smaller values of k. The rows of the matrix are then stably deduplicated. This removes two types of duplicates: 1. Pairs that have been found more than once for the current k (such as 3,4, which also results from 6,8 when dividing by its gcd); 2. Pairs that were already found with smaller k. In fact, each iteration k finds all pairs that were already found for previous iterations. But it may find them in a different order. For example, k=25 will find the triple 7,24,25 and not 20,21,29 (because c cannot exceed k). Later on, iteration k=29 will find both, but with 20,21,29 before 7,24,25 (order is increasing b, then a). This is why, instead of keeping all pairs found for the latest k, we append them to the previous ones and deduplicate stably. Doing so assures that the order is the same for any input n. The above guarantees that each primitive Pythagorean triple will eventually appear, and it will appear only once. For input n, the loop over k finishes when at least n valid triples have been obtained; and then the n-th triple is ouput. Try it online! Or use this modified code to see the first n triples: @:Ut&+wmR&fZd1Mhw/vXutnGE<]G:Y)tU&2sX^h  Try it online! • Nice explanation. – AdmBorkBork Jul 25 '17 at 15:47 # Haskell, 98 bytes f 0=(3,4,5) f n|let(a,b,c)=f$div(n-1)3;r=mod n 3;d=a*(-1)^2^r;e=b*(-1)^r;s=2*(d+e+c)=(s-d,s-e,s+c)


Try it online!

### How it works

This interprets the bijective base-3 digits of n as a path down the tree of primitive Pythagorean triples. It runs with no searching in O(log n) operations. # Jelly, 19 18 bytes

*g/²_/
4*œc3UṢÇÐḟḣ


The program takes a 1-based index n and prints the first n PPTs [c, b, a] in lexicographical order.

This is a O(64n) solution, so TIO will choke on inputs 4 and higher. I'll work on making it faster.

Try it online!

### Alternate version, O(n3), probably valid

To find the nth triplet – [cn, bn, an] – the solution above assumes that cn ≤ 4n, which is easy to verify. However, A020882 proves that cn ~ 2πn, so there is a k such that cn ≤ kn for all n.

If we can take k = 7, the solution below is also valid (and much, much faster).

*g/²_/
×7œc3UṢÇÐḟḣ


Try it online!

### How it works

4*œc3UṢÇÐḟḣ  Main link. Argument: n

4*           Compute 4**n, the n-th power of 4.
œc3        Take all 3-combinations of the set {1, ..., 4**n}, each sorted in
ascending order. The triples themselves are sorted lexicographically.
U       Upend; reverse each triple [a, b, c], yielding [c, b, a].
Ṣ      Sort the descending triples lexicographically. This ensures that their
order is independent of n.
ÇÐḟ   Keep only triples for which the helper link returns a falsy value.
ḣ  Dyadic head; take the first n triples.

*g/²_/       Helper link. Argument: [c, b, a]

g/          Reduce [c, b, a] by greatest common divisor, yielding g.
*            Elevate the integers to that power, computing [c**g, b**g, a**g].
²         Square, yielding [c**2g, b**2g, a**2g].
_/       Reduce by subtraction, yielding c**2g - b**2g - a**2g.
Fermat's Last Theorem assures that this difference will be non-zero
whenever g > 1, so this yields 0 iff g = 1 and c² = a² = b².


## JavaScript (ES7), 106105 103 bytes

Outputs the Nth PPT. Results are 1-indexed and ordered by the value of b.

n=>(g=(a,b)=>b?g(b,a%b):a,F=a=>(x=a*a+b*b,c=x**.5|0)*c-x*g(a,g(b,c))||--n?F(a-b?a+1:!++b):[a,b,c])(b=1)


### Demo

let f =

n=>(g=(a,b)=>b?g(b,a%b):a,F=a=>(x=a*a+b*b,c=x**.5|0)*c-x*g(a,g(b,c))||--n?F(a-b?a+1:!++b):[a,b,c])(b=1)

for(n = 1; n <= 20; n++) {
console.log('#' + n + ' : ' + f(n));
}

# MATL, 63 bytes

3lvi:"t"[HlHllO;aOlOHl]!@Y*2eh]]!XuGY)&*tt[lO;Oa]*ssD2)ED2Xy*ss


Try it online!

A lesson in golfing gone horribly wrong. I'm posting it anyway because I'm wondering if there are ways to golf this better.

I based this on this Wikipedia page in combination with Euclid's formula, to constructively generate all triples rather than trial-and-error approaches.

First, odd coprime pairs are generated as a ternary tree. This is done as a big matrix multiplication, accounting for most of the byte count. Then, Euclid's formula is applied, possibly also in a very byte-wasting manner. If anyone has tips for these two parts, I'd love to hear them.

Note that, to save bytes, this program generates a tree the same depth as the input, rather than log3(n). Also, children are generated for each row rather than only for the last row of the tree, and then filtered again with Xu. So much for an efficient constructive approach.

3lv % Push root node of ternary tree
i:" % Generate a tree of depth of input (WAY too large, but golfy)
t"  % loop over all nodes (golfier than looping over last tree row)
[HlHllO;aOlOHl]! % Matrix to generate three children of current node
@Y* % Multiply with current node to get children
2e  % Reshape to get node pairs
h]] % Append to tree, exit loops
!Xu % Remove duplicates (more efficient to do it before last ] but golfier this way)
GY) % Select n-th odd coprime pair
&*tt % Multiply with it's own transpose to get [m²,m*n;m*n,n²]
[lO;Oa]*ssD % Sum of matrix multiplication = m²-n² to get a
2)ED % Second element doubled for b=2mn
2Xy*ss % Sum of matrix multiplication with identify matrix to get c=m²+n²


## Haskell, 65 bytes

([(a,b,c)|c<-[5..],b<-[1..c],gcd c b<2,a<-[1..b],a^2+b^2==c^2]!!)


0-based indexing. For a given c, b runs up to c and a up to b, so c > b > a always holds.

Try it online!

# Python, 675048 46 Bytes

Using the formulas found on wikipedia,

a=m*n, b=(m^2-n^2)/2, c=(m^2+n^2)/2

where m>n>0 and m and n are coprimes and odd. Here is the code

lambda n:[3+2*n,~-(3+2*n)**2-1/2,-~(3+2*n)**2/2]


-17 bytes thanks to @Martin Ender

Try it Online

Works by always having the value of the n variable in the equation being 1, meaning that m simply is any other odd value, in this case, 3+2*n where n is the number of the primitive Pythagorean triple. This allows us to assume the value of 1 for all n values.

• Welcome to PPCG! Unnamed functions are fine, so you don't need to assign the lambda to a (and if you did, you could get rid of the two spaces there). I'm also not sure why you print there, you could just return the values from the lambda itself. – Martin Ender Jul 26 '17 at 19:50
• "you can prove (informally is fine) that your algorithm can generate every possible unique PPT". But this method only generates those where the hypotenuse is 1 longer than a leg. It never generates 8,15,17, for example. – Rosie F Aug 1 at 6:33

# Husk, 18 bytes

↑üOf§=F⌋ȯ¬Ḟ-m□ΠR3N


Try it online!

-4 bytes thank to Zgarb, with inspiration from Dennis

Super-slow brute force approach, won't work on TIO for inputs larger than 1. You can try a more manageable version, limited to a,b≤200 here

### Explanation

↑üOf§=F⌋ȯ¬Ḟ-m□ΠR3N
ΠR3N   Get all triples of natural numbers
f                 Keep only those triples where
F⌋                their GCD
§=                  is equal to
ȯ¬Ḟ-m□          the logical negation of c²-b²-a²
üO                  Remove duplicates by their sorted version
↑                    Get the first <input> values of this sequence

• 20 bytes by combining the map and filter, even slower. – Zgarb Jul 26 '17 at 9:19
• @Zgarb thank you! I managed to golf an extra byte :) – Leo Jul 26 '17 at 11:29
• 18 bytes with the subtraction trick from Dennis's Jelly answer. – Zgarb Jul 26 '17 at 12:59
• @Zgarb nice! Though I'm having a doubt: could there be two different triples with the same c? in that case this solution would need to be fixed – Leo Jul 26 '17 at 17:42
• Hmm, actually there are many triples with the same c. This 18-byte solution (which uses another trick of Dennis) works regardless. – Zgarb Jul 26 '17 at 18:21

# Mathematica, 89 bytes

using Solve ordered by c

SortBy[{a,b,c}/.Solve[a^2+b^2==c^2&&GCD[a,b]==1&&0<a<b<c<9#,{a,b,c},Integers],Last][[#]]&


# Mathematica, 124 bytes

(s={};Table[If[IntegerQ[c=Sqrt[a^2+b^2]]&&GCD[a,b]==1,AppendTo[s,{a,b,c}]],{a,9#},{b,9#}];SortBy[Union[Sort/@s],Last][[#]])&


# R (+ numbers), 88 bytes

n=scan();while(all(nrow(T)<n))T=numbers::pythagorean_triples(5,5+(F<-F+1));T[n,3:5]


For using a builtin to generate the numbers, it takes actually a surprising amount of bytes to get what we want. The builtin takes two arguments c1 and c2, and returns those triplets that have c >= c1 & c <= c2. This makes it slightly annoying to get to the n-th triplet. This will just keep increasing c2 1 at a time until the output is just enough rows.

# PHP, 273 Bytes

function t($n){$x=[];for($c=3;;$c++)for($b=2;$b<$c;$b++)for($a=2;$a<$b;$a++)if(d($a,$b,$c)&&$a**2+$b**2==$c**2){$x[]=[$a,$b,$c];if(--$n==0)return$x;}}function d($a,$b,$c){for($i=2;$i<$a;$i++)if($a%$i==0&&$b%$i==0||$a%$i==0&&$c%$i==0||$b%$i==0&&$c%$i==0)return 0;return 1;}  • t($n) returns an array of [a,b,c] with ordering a < b < c
• Returns a zero-based index

Try it online! (code is readable there too)

# C, 158 bytes

I believe this is my first submission here, so you can most probably do better.

#include<stdio.h>
void f(n){int i=0,j=3,k,l;while(1){for(k=1;k<j;k++){for(l=k;l<j;l++){if(j*j==k*k+l*l)i++;if(i==n){printf("%d %d %d",j,k,l);return;}}}j++;};}


And ungolfed version:

#include <stdio.h>

void f(n)
{
int i=0, j=3, k, l;
while (1) {
for (k=1; k<j; k++) {
for (l=k; l<j; l++) {
if (j*j==k*k+l*l)
i++;
if (i==n) {
printf("%d %d %d\n", j, k, l);
return;
}
}
}
j++;
};
}

void main()
{
int i;

scanf("%d", &i);

f(i);
printf("\n");
}


For a2 + b2 = c2, the ordering is increasing c then increasing a.

There can no be twice the same PPT as b is at leas a in this algorithm.

• Welcome to PPCG! – JAD Jul 26 '17 at 11:31

# Jelly, 27 25 bytes

⁽0(ḃs
Ɠḃd2Ḥ’×€Ç
3r5DṭÇæ×/


This is an implementation of the tree approach from @AndersKaseorg's Haskell answer, with a different branch order. The program uses 0-based indexing and takes input from STDIN.

Try it online!

### Background

As mentioned on the Wikipedia page Tree of primitive Pythagorean triples, every PPT can be obtained by repeatedly left-multiplying the row vector (3, 4, 5) by matrices with certain properties.

In each iteration, the previous result can be left-multiplied by either A, B, or C, which can be chosen as follows. When A, B, and C are fixed, each PPT can be obtained in a unique fashion.

### How it works

3r5DṭÇæ×/  Main link. No arguments.

3          Set the argument and the return value to 3.
r5        Create a range from 3 to 5, i.e., [3, 4, 5].
D       Decimal; convert each integer to base 10, yielding [, , ].
Ç     Call the second helper link with argument 3.
ṭ      Tack; append [, , ] to the result.
æ×/  Reduce by matrix multiplication.

Ɠḃd2Ḥ’×€Ç  Second helper link. Argument: 3

Ɠ          Read and evaluate one line of input, yielding an integer n.
ḃ         Convert n to bijective base 3.
d2       Divmod 2; map each digit d to [d/2, d%2].
Ḥ      Unhalve; multiply the results by 2.
’     Decrement the doubled results.
The previous four atoms apply the following mapping to the digits.
1 -> [0, 1] -> [0, 2] -> [-1,  1]
2 -> [1, 0] -> [2, 0] -> [ 1, -1]
3 -> [1, 1] -> [2, 2] -> [ 1,  1]
Ç  Call the helper link with argument 3, yielding the following 2D array.
[[ 1,  2,  2],
[ 2,  1,  2],
[ 2,  2,  3]]
×€   Multiply each [-1,  1], [ 1, -1], and [ 1,  1] by that matrix, using
vectorizing multiplication (not matrix multiplication), yielding one
of the following three 2D arrays.

[[-1,  2,  2],    [[ 1, -2,  2],    [[ 1,  2,  2],
[-2,  1,  2],     [ 2, -1,  2],     [ 2,  1,  2],
[-2,  2,  3]]     [ 2, -2,  3]]     [ 2,  2,  3]]

⁽0(ḃs      First helper link. Argument: 3

⁽0(        Numeric literal; yield 13041.
ḃ       Convert 13041 to bijective base 3, yielding [1, 2, 2, 2, 1, 2, 2, 2, 3].
s      Split the result into chunks of length 3, yielding the aforementioned
2D array.


# APL(NARS), 90 chars, 180 bytes

{a⊃⍨⍵⊃⍋↑¨a←{⍵[⍋⍵]}¨a/⍨{1=∨/⍵}¨a←{(-/k),(×/2,⍵),+/k←⍵*2}¨a/⍨{>/⍵}¨a←,a∘.,a←⍳(⌊2⍟2+⍵)×9+⌊√⍵}


if the argument of the function above is ⍵, the above function would return the element of the index ⍵ (1 based) of the array has elements Pythagorean triples (a,b,c) where a<=b<=c and that array is order first for a, (the side more short), then for b (the other side not hypotenuse). There would be something wrong because there is not seen where I order for b too... test:

  f←{a⊃⍨⍵⊃⍋↑¨a←{⍵[⍋⍵]}¨a/⍨{1=∨/⍵}¨a←{(-/k),(×/2,⍵),+/k←⍵*2}¨a/⍨{>/⍵}¨a←,a∘.,a←⍳(⌊2⍟2+⍵)×9+⌊√⍵}
f¨1..10
3 4 5  5 12 13  7 24 25  8 15 17  9 40 41  11 60 61  12 35 37  13 84 85  15 112 113  16 63 65


it is related with http://oeis.org/A020884 and http://oeis.org/A020884/b020884.txt

A020884: Ordered short legs of primitive Pythagorean triangles.

  ↑¨f¨1..23
3 5 7 8 9 11 12 13 15 16 17 19 20 20 21 23 24 25 27 28 28 29 31
f 999
716 128163 128165
f 1000
717 28556 28565


I don't know if it is right, it seems function give correct result of first side of triangle until 1000 but I don't know for the remain, and possible could be some triple not right too even <1000.

# JavaScript, 101 bytes

By Euclid's formula all primitive Pythagorean triples can be generated from integers m and n with m>n>0, m+n odd gcd(m,n)==1 (Wikipedia)

This function enumerates all m,n pair incrementing m starting from m=2 and decrementing n by 2 starting from m-1 (so that m+n is odd)

c=>eval("g=(a,b)=>b?g(b,a%b):a;for(m=2,n=1;c-=g(m,n)<2;(n-=2)>0||(n=m++));[m*m-n*n,2*m*n,m*m+n*n]")


Less golfed

c => {
g = (a,b) => b ? g(b,a%b) : a;
for( m = 2, n = 1;
g(m,n) < 2 ? --c : c;
(n -= 2) > 0 || (n = m++))
/* empty for body */;
return [m*m - n*n, 2*m*n, m*m + n*n]
}


Test

F=
c=>eval("g=(a,b)=>b?g(b,a%b):a;for(m=2,n=1;c-=g(m,n)<2;(n-=2)>0||(n=m++));[m*m-n*n,2*m*n,m*m+n*n]")

for(i=1;i<=50;i++) console.log(i+' '+F(i))