%I #81 Jan 03 2023 21:06:07
%S 3,4,5,21,20,29,105,88,137,465,368,593,1953,1504,2465,8001,6080,10049,
%T 32385,24448,40577,130305,98048,163073,522753,392704,653825,2094081,
%U 1571840,2618369,8382465,6289408,10479617,33542145,25161728,41930753,134193153,100655104,167747585,536821761,402636800,671039489,2147385345,1610579968,2684256257
%N Three-column array pPT read by rows: subsequence of primitive Pythagorean triples (x, y, z) with x = A153893^2 - A000079^2, y = 2*A153893*A000079, z = A153893^2 + A000079^2, ordered by increasing z.
%C Let [h21] = {{1, 3}, {0, 2}} be the matrix [h_2]*[h_1] in Firstov's notation, from eqs. (24) and (39). Then primitive Pythagorean triples (pPT) (x(n), y(n), z(n)) = (u(n)^2 - v(n)^2, 2*u(n)*v(n), u(n)^2 + v(n)^2), with u(n) and v(n) of different parity, gcd(u(n), v(n)) = 1, and u(n) > v(n) > 0, are generated by (u(n), v(n))^T = [h21]^n*(2,1)^T (T for transpose).
%C For n > 0: (x(n), y(n), z(n)) = (1, 0, 1) (mod 4). Thus some z are Pythagorean primes (A002144).
%C The triples converge to the proportion (4:3:5) with:
%C lim_{n->infinity} x(n)/y(n) = 4/3, lim_{n->infinity} y(n)/z(n) = 3/5.
%C Altitude h(n) = x(n)*y(n)/z(n) is an irreducible fraction because of primitivity.
%C From _Wolfdieter Lang_, Jun 13 2020: (Start)
%C [h21]^n = sqrt(2)^n*(S(n, 3/sqrt(2))*[1_3] + S(n-1, 3/sqrt(2))*(1/sqrt(2))*([h21] - 3*[1_3])) with the Chebyshev S polynomials (A049310).
%C u(n) = sqrt(2)^n*(2*S(n, 3/sqrt(2)) - (1/sqrt(2))*S(n-1, 3/sqrt(2)))
%C = A153893(n),
%C v(n) = sqrt(2)^n*(S(n, 3/sqrt(2)) - (1/sqrt(2))*S(n-1, 3/sqrt(2)))
%C = A000079(n). Proof from the recurrence, using the Cayley-Hamilton theorem.
%C With the monic Chebyshev T polynomials, called R in A127672:
%C x(n)/3 = 2^(n+1)*(R(2*(n+1), 3/sqrt(2)) - (sqrt(2)/3)*R(2*n+1,3/sqrt(2)) - 1) = A171477(n),
%C y(n)/4 = 3*2^(n-1)*(sqrt(2)*R(2*n+1,3/sqrt(2)) - R(2*n,3/sqrt(2)) - 1/3)
%C = A010036(n),
%C z(n) = 3*2^(n+1)*((3/sqrt(2))*R(2*n+1, 3/sqrt(2)) - (4/3)*R(2*n,3/sqrt(2)) - 1).
%C Using 2^n*Rnx(2*n, 3/sqrt(2)) = A052539(n) = 2^(2*n) + 1, and
%C 2^(n)*(sqrt(2)/3)*Rnx(2*n+1, 3/sqrt(2)) = A007583(n) = (2^(2*n + 1) + 1)/3,
%C produces the explicit formulas given by the author in the formula section.
%C G.f.s for {x(n)} G0(x) = 3/((1 - 4*x)*(1 - 2*x)*(1 - x)), for {y(n)} G1(x) = 4*(1-x)/((1 - 4*x)*(1 - 2*x)), and for {z(n)} = (5 - 6*x + 4*x^2)/((1 - 4*x)*(1 - 2*x)*(1 - x)). This produces the g.f. for the array, read as sequence {a(n)}: G(x) = G0(x^3) + x*G1(x^3) + x^2*G2(x^3) given in the formula section by _Colin Barker_.
%C (End)
%H V. E. Firstov, <a href="http://mi.mathnet.ru/eng/mz4074">A Special Matrix Transformation Semigroup of Primitive Pairs and the Genealogy of Pythagorean Triples</a>; Mathematical Notes, volume 84, number 2, August 2008, pages 263-279; Link of the page (for the Russian article).
%H Ralf Steiner, <a href="https://www.researchgate.net/publication/341131567_Weitere_spezielle_Folge_primitiver_pythagoraischer_Dreiecke">Weitere spezielle Folge primitiver pythagoraischer Dreiecke</a>, ResearchGate.
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples">Tree of primitive Pythagorean triples</a>.
%H <a href="/index/Rec#order_09">Index entries for linear recurrences with constant coefficients</a>, signature (0,0,7,0,0,-14,0,0,8).
%F The three-column array PT(n, k) is for k = 0, 1, 2: x(n), y(n), z(n), for n >= 0, with
%F x(n) = a(3*n + 0) = A153893(n)^2 - A000079(n)^2 = 1 - 3*2^(n+1) + 2^(2*n+3) = binomial(2^(n+2) - 1, 2) = 3*A171477(n),
%F y(n) = a(3*n + 1) = 2*A153893(n)*A000079(n) = 2^(n+1)*(-1 + 3*2^n) = 4*A010036(n),
%F z(n) = a(3*n + 2) = A153893(n)^2 + A000079(n)^2 = 1 - 6*2^n + 10*2^(2*n).
%F From _Colin Barker_, May 08 2020: (Start)
%F G.f. (read as sequence {a(n)}): (3 + 4*x + 5*x^2 - 8*x^4 - 6*x^5 + 4*x^7 + 4*x^8) / ((1 - x)*(1 + x + x^2)*(1 - 2*x^3)*(1 - 4*x^3)).
%F a(n) = 7*a(n-3) - 14*a(n-6) + 8*a(n-9), for n > 8.
%F (End)
%e The three-column array pPT(n,k) begins:
%e n\k 0 1 2
%e -------------------------------
%e 0: 3 4 5
%e 1: 21 20 29
%e 2: 105 88 137
%e 3: 465 368 593
%e 4: 1953 1504 2465
%e 5: 8001 6080 10049
%e 6: 32385 24448 40577
%e 7: 130305 98048 163073
%e 8: 522753 392704 653825
%e 9: 2094081 1571840 2618369
%e 10: 8382465 6289408 10479617
%e ... - _Wolfdieter Lang_, Jun 13 2020
%t h21={{1, 3}, {0, 2}}; l = {}; Do[v = MatrixPower[h21, n, {2, 1}]; p = v[[1]]; q = v[[2]];
%t a = p^2 - q^2; b = 2 p q; c = p^2 + q^2; l = AppendTo[l, {a, b, c}], {n, 0, 14}]; l // Flatten
%o (PARI) Vec((3 + 4*x + 5*x^2 - 8*x^4 - 6*x^5 + 4*x^7 + 4*x^8) / ((1 - x)*(1 + x + x^2)*(1 - 2*x^3)*(1 - 4*x^3)) + O(x^35)) \\ _Colin Barker_, Jun 12 2020
%Y Cf. A000079, A010035, A049310, A083329, A093357, A010036, A134057, A153893, A171477.
%K nonn,easy,tabf
%O 0,1
%A _Ralf Steiner_, May 07 2020
%E Edited, and corrected proportion by _Wolfdieter Lang_, Jun 13 2020
%E Minor grammatical edits. - _N. J. A. Sloane_, Sep 12 2020