Pythagorean Triple Sequence

Question

A Pythagorean triple consists of three positive integers a, b, and c, such that a² + b² = c². 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)

Answer 1

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 {5², 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 {5², 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

Answer 2

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

Answer 3

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².

Answer 4

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!

Answer 5

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)

Answer 6

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

Answer 7

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)

Answer 8

J, 54 47 bytes

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

TIO

greedy split of prime factors into coprime factors

old 54 bytes TIO

Answer 9

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!

Answer 10

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}

Answer 11

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++);

Answer 12

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)². 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.

Answer 13

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

Answer 14

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