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!):


A018928
A077034 (the odd terms)
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