# Pythagorean Triple Sequence

A Pythagorean triple consists of three positive integers a, b, and c, such that a2 + b2 = c2. Such a triple is commonly written (a, b, c), and a well-known example is (3, 4, 5). If (a, b, c) is a Pythagorean triple, then so is (ka, kb, kc) for any positive integer k. A primitive Pythagorean triple is one in which a, b and c are coprime.

Using this knowledge, we can create a sequence by chaining together the least lengths of triples, where the next element in the sequence is the hypotenuse (largest number) of the smallest primitive Pythagorean triple containing the previous element as the smallest one of its lengths.

Start with the smallest primitive Pythagorean triple (3, 4, 5). The sequence begins with 3, and the hypotenuse (next element in the sequence) is 5. Then find the smallest primitive Pythagorean triple with 5 as a leg, and you get (5, 12, 13). So the sequence continues with 13.

Either output the sequence forever, or take an integer input n and output the first n elements of the sequence, either zero or one indexed.

You need to support output at least through and including 28455997, but if the limit of the data type you are using was suddenly raised, it would need to work for that new limit. So you cannot hard code a list of numbers.

3
5
13
85
157
12325
90733
2449525
28455997
295742792965
171480834409967437
656310093705697045
1616599508725767821225590944157
4461691012090851100342993272805
115366949386695884000892071602798585632943213
12002377162350258332845595301471273220420939451301220405


OEIS A239381

Similar sequences (don't output these!):

• Is there any time limit? – Loovjo Oct 28 '16 at 15:59
• @Loovjo No, but you should know/prove that your output is correct. There are some similar sequences where the output differs after 12325. – mbomb007 Oct 28 '16 at 16:01
• The similar sequence I'm thinking of differs after 85... its next term is 3613 (can you guess what it is yet?) – Neil Oct 28 '16 at 16:07
• @Neil A quick search yields A053630, the Pythagorean spiral. I referenced the two in the challenge though, because while working on making my implementation I accidentally reached those two sequences or similar to them. – mbomb007 Oct 28 '16 at 16:18
• Indeed, had I been more awake I could have just looked it up myself... – Neil Oct 28 '16 at 16:27

# Jelly, 19 bytes

o3ṄÆF*/€ŒPP€²+Ṛ$HṂß  Saved a byte thanks to @Dennis by refactoring to an infinite sequence. Takes no input and arguments, then outputs the sequence infinitely by printing each term as it computes them. This method slows down as the numbers get larger since it depends on prime factorization. Try it online! This calculates the next term by computing the prime power factorization of the current term. For 12325, this is {52, 17, 29}. There is a variant of Euclid's formula for calculating Pythagorean triples {a, b, c}, where m > n and the triple is primitive iff m and n are coprime. To calculate the next primitive root from 12325, find m and n such that mn = 12325 and choose m,n so that gcd(m, n) = 1. Then generate all pairs of m,n by creating all subsets of {52, 17, 29} and finding the product of each of those subsets which are {1, 25, 17, 29, 425, 725, 493, 12325}. Then divide 12325 by each value and pair so that each pair is m,n. Compute the formula for c using each pair and take the minimum which is 90733. • The previous method failed for determining the next term after 228034970321525477033478437478475683098735674620405573717049066152557390539189785244849203205. The previous method chose the last value as a factor when the correct choice was the 3rd and last primes. The new method is slower but will always work since it tests all pairs of coprimes to find the minimal hypotenuse. ## Explanation o3ṄÆF*/€ŒPP€²+Ṛ$HṂß  Main link. Input: 0 if none, else an integer P
o3                   Logical OR with 3, returns P if non-zero else 3
Ṅ                  Println and pass the value
ÆF                Factor into [prime, exponent] pairs
*/€             Reduce each pair using exponentation to get the prime powers
ŒP           Powerset of those
P€         Product of each
²        Square each
$Monadic chain + Add vectorized with Ṛ the reverse H Halve Ṃ Minimum ß Call recursively on this value  • Wow, this is really fast! – mbomb007 Oct 28 '16 at 19:12 • o3ṄÆfµṪ,P²SHß with infinite output saves a byte. – Dennis Oct 28 '16 at 19:29 # Brachylog, 36 bytes 3{@wB:?>:^a+~^=C:B:?:{$pd}ac#d,C:1&}


Try it online!

You have to wait for the program to time out (1 minute) before TIO flushes the output. In SWI-Prolog's REPL this prints as soon as it finds the value.

This will print the sequence forever.

After a few minutes on the offline SWI-Prolog's interpreter, I obtained 90733 after 12325. I stopped it after this point.

This isn't complete bruteforce as it uses constraints to find pythagorean triples, though it is obviously not optimised for speed.

### Explanation

3{                                 }    Call this predicate with 3 as Input
@w                                    Write the Input followed by a line break
B:?>                                B > Input
+                            The sum...
:^a                             ...of Input^2 with B^2...
~^                          ...must equal a number which is itself a square
=C                        Assign a fitting value to that number and call it C
C:B:?:{$pd}a Get the lists of prime factors of C, B and Input without duplicates c#d, Concatenate into a single list; all values must be different C:1& Call recursively with C as Input  # Perl, 73 bytes for($_=3;$_<1e9;$_=$a**2+$b**2){$a++until($b=($_+$a**2)**.5)==($b|0);say}  All Pythagorean Triples a²+b²=c² satisfy a=r(m²-n²), b=2rmn, c=r(m²+n²) for some integers r,m,n. When r=1 and m,n are coprime with exactly one being divisible by 2, then a,b,c is a primitive triple, where a,b,c are all pairwise coprime. With this in mind, given some a, I use a brute-force algorithm to calculate the smallest n such that a²-n² is a square, namely m². Then, c is equal to n²+m². • Possible typo in your explanation: you search for n such that a+n² is a square. – Neil Oct 29 '16 at 15:49 # Python 3, 178 bytes from math import* p,n=[3,5],int(input()) while len(p)<n: for i in range(p[-1],p[-1]**2): v=sqrt(i**2+p[-1]**2) if v==int(v)and gcd(i,p[-1])==1: p+=[int(v)];break print(p)  This is basically just a brute force algorithm, and thus is very slow. It takes the amount of terms to output as the input. I'm not 100% sure about the correctness of this algorithm, the program checks for the other leg up to the first leg squared, which I believe is enough, but I haven't done the math. Try it on repl.it! (Outdated) (Please don't try it for numbers greater than 10, it will be very slow) • You can switch to Python 3.5 and use math.gcd. Also, use p+=[...] instead of p.append(...). And <2 instead of ==1. And the if can all be on one line. – mbomb007 Oct 28 '16 at 16:38 • You can still do the last 2 improvements I suggested. – mbomb007 Oct 28 '16 at 18:14 • 161 bytes – Mr. Xcoder Jul 25 '17 at 13:54 • Loovjo, you gonna golf your code using the suggestions or not? – mbomb007 May 13 '18 at 3:14 # MATL, 27 bytes Ii:"I@Yyt1\~?3MZdZdq]}6MXI  This produces the first terms of the sequence. Input is 0-based. The code is very inefficient. The online compiler times out for inputs greater than 5. Input 6 took a minute and a half offline (and produced the correct 90733 as 6-th term). Try it online! I % Push 3 (predefined value of clipboard I) i % Input n :" % For each (i.e. execute n times)  % Do...while I % Push clipboard I. This is the latest term of the sequence @ % Push iteration index, starting at 1 Yy % Hypotenuse of those two values t1\ % Duplicate. Decimal part ~? % If it is zero: we may have found the next term. But we still % need to test for co-primality 3M % Push the two inputs of the latest call to the hypotenuse % function. The stack now contains the hypotenuse and the % two legs ZdZd % Call GCD twice, to obtain the GCD of those three numbers q % Subtract 1. If the three numbers were co-prime this gives % 0, so the do...while loop will be exited (but the "finally" % part will be executed first). If they were not co-prime % this gives non-zero, so the do...while loop proceeds % with the next iteration ] % End if % If the decimal part was non-zero: the duplicate of the % hypotenuse that is now on the top of the stack will be used % as the (do...while) loop condition. Since it is non-zero, % the loop will proceed with the next iteration } % Finally (i.e. execute before exiting the do...while loop) 6M % Push the second input to the hypotenuse function, which is % the new term of the sequence XI % Copy this new term into clipboard I % Implicitly end do...while % Implicitly end for each % Implicitly display the stack, containing the sequence terms  ## Racket 106 bytes (let p((h 3))(println h)(let p2((i 1))(define g(sqrt(+(* h h)(* i i))))(if(integer? g)(p g)(p2(add1 i)))))  Ungolfed: (define (f) (let loop ((h 3)) (let loop2 ((i 1)) (define g (sqrt (+(* h h) (* i i)))) (if (not(integer? g)) (loop2 (add1 i)) (begin (printf "~a ~a ~a~n" h i g) (loop g))))))  Testing: (f)  Output of golfed version: 3 5 13 85 157 12325 12461 106285 276341 339709 10363909 17238541  Output of ungolfed version: 3 4 5 5 12 13 13 84 85 85 132 157 157 12324 12325 12325 1836 12461 12461 105552 106285 106285 255084 276341 276341 197580 339709 339709 10358340 10363909 10363909 13775220 17238541  (Error after this on my machine) • The golfed code prints out only hypotenuse of sequence. Ungolfed versions shows all three to clarify triplets not mentioned in question. – rnso Oct 30 '16 at 4:14 # Wolfram Language (Mathematica), 74 bytes Nest[c/.#&@@Solve[{#^2+b^2==c^2,#<b<c&&GCD[b,c]==1},{b,c},Integers]&,3,#]&  Try it online! # Wolfram Language (Mathematica), 74 bytes Nest[#&@@Cases[Divisors@#,d_/;GCD[d,#/d]<2&&d^2>#:>(d^2+(#/d)^2)/2]&,3,#]&  Try it online! # PHP, 139 bytes for($k=3;$i=$k,print("$k\n");)for($j=$i+1;($k=sqrt($m=$i*$i+$j*$j))>(int)$k||gmp_intval(gmp_gcd(gmp_gcd((int)$i,(int)$j),(int)$k))>1;$j++);


The above code breaks after 28455997 on 32-bit systems. If higher numbers are needed, it becomes 156 bytes:

for($k=3;$i=$k,print("$k\n");)for($j=$i+1;!gmp_perfect_square($m=bcadd(bcpow($i,2),bcpow($j,2)))||gmp_intval(gmp_gcd(gmp_gcd($i,$j),$k=bcsqrt($m)))>1;$j++);


# Java 8, 133 Bytes

-25 bytes thanks to miles Using n*n instead of Math.pow(n, 2)

-24 bytes thanks to miles Using for loops instead of while, changing datatype, eliminating () due to order of operations

()->{long b=3,c,n;for(;;){for(n=1;;n++){c=b+2*n*n;double d=Math.sqrt(c*c-b*b);if(d==(int)d&b<d){System.out.println(b);break;}}b=c;}};


Uses the fact that

for any pair of integers m > n > 0. Therefore, C is equal to A plus 2(N)2. The function above finds the least value of N that satisfies this relation, while making the second element of the Pythagorean triple an integer and greater than the first element. Then it sets the value of the first element to the third element and repeats with the updated first element.

Ungolfed:

void printPythagoreanTriples() {
long firstElement = 3, thirdElement, n;
while (true) {
for (n = 1; ; n++) {
thirdElement = firstElement + (2 * n * n);
double secondElement = Math.sqrt(thirdElement * thirdElement - firstElement * firstElement);
if (secondElement == (int) secondElement && firstElement < secondElement) {
System.out.println("Found Pythagorean Triple [" +
firstElement + ", " +
secondElement + ", " +
thirdElement + "]");
break;
}
}
firstElement = thirdElement;
}
}


Ideone it!

*The ideone does not print the last required element due to time limits, however as you can see through the logic of the program and the ungolfed version (which prints the 28455997 as the third element of the previous Pythagorean triple rather than the first element of the next), the values are, with a higher time limit, printed.

• Couldn't you use n*n instead of Math.pow(n,2)? – miles Oct 29 '16 at 22:55
• I don't know why I didn't think of that... I'll add that right away. Thank you @miles – Mario Ishac Oct 29 '16 at 22:58
• I shaved some more off using for loops to get it down to 133 bytes ()->{long b=3,c,n;for(;;){for(n=1;;n++){c=b+2*n*n;double d=Math.sqrt(c*c-b*b);if(d==(int)d&b<d){System.out.println(b);break;}}b=c;}}; – miles Oct 29 '16 at 23:33

# Python 3.5, 97 bytes

Wrong output after 28455997, because of the limits of the floating point data type. The sqrt function isn't good enough, but if the precision was magically increased, it'd work.

Pretty simple to understand. Incrementing c by two instead of one cuts the runtime in half, and only odd numbers need to be checked anyway, because the elements are always odd.

import math
c=a=3
while 1:
c+=2;b=(c*c-a*a)**.5;i=int(b)
if math.gcd(a,i)<2<a<b==i:print(a);a=c


Try it online

The program cannot be run on Ideone, because Ideone uses Python 3.4

For output to stay accurate longer, I'd have to use decimal:

import math
from decimal import*
c=a=3
while 1:
c+=2;b=Decimal(c*c-a*a).sqrt();i=int(b)
if i==b>a>2>math.gcd(a,i):print(a);a=c


Try it online

To stay accurate indefinitely, I could do something horrid like this (increasing the precision required every single iteration:

import math
from decimal import*
c=a=3
while 1:
c+=2;b=Decimal(c*c-a*a).sqrt();i=int(b);getcontext().prec+=1
if i==b>a>2>math.gcd(a,i):print(a);a=c


# J, 54 47 bytes

-:@+/@:*:@((*/:~)/)@,.&1@x:@(^/)@(2&p:)^:(<12)3


TIO

greedy split of prime factors into coprime factors

old 54 bytes TIO

# Pari/GP, 71 bytes

n->a=3;for(i=1,n,a=[d^2+e^2|d<-divisors(a),gcd(d,e=a/d)<2&&d>e][1]/2);a


Try it online!

# APL(NARS), 169 chars, 338 bytes

h←{{(m n)←⍵⋄(mm nn)←⍵*2⋄(2÷⍨nn+mm),(2÷⍨nn-mm),m×n}a⊃⍨b⍳⌊/b←{⍵[2]}¨a←a/⍨{(≤/⍵)∧1=∨/⍵}¨a←(w÷a),¨a←∪×/¨{k←≢b←1,π⍵⋄∪{b[⍵]}¨↑∪/101 1‼k k}w←⍵}⋄p←{⍺=1:⍵⋄⍵,(⍺-1)∇↑h ⍵}⋄q←{⍵p 3x}


test ok until 14 as argument of q:

  q 1
3
q 2
3 5
q 10
3 5 13 85 157 12325 90733 2449525 28455997 295742792965
q 12
3 5 13 85 157 12325 90733 2449525 28455997 295742792965 171480834409967437 656310093705697045
q 13
3 5 13 85 157 12325 90733 2449525 28455997 295742792965 171480834409967437 656310093705697045
1616599508725767821225590944157
q 14
NONCE ERROR
q 14
∧


this below would find all divisors of its argument...

∪×/¨{k←≢b←1,π⍵⋄∪{b[⍵]}¨↑∪/101 1‼k k}


# JavaScript (Node.js), 101 bytes

_=>{for(var a=3,b,c;;){for(c=1;;c++)if(b=a+2*c*c,d=Math.sqrt(b*b-a*a),d==+d&a<d){alert(a);break}a=b}}


Try it online!

Suggestions on golfing are welcome