Starting with $(165,52,173)$, we attempt to find the $p,q$ pair which generates this triple, i.e., $p^2-q^2=165,2pq=52,p^2+q^2=173$. Clearly, we have $2q^2=173-165=8\implies q=2$ and thus $p=13$.
The ancestor of this triple arises from $(|p-2q|,q)\text{ or }(q,|p-2q|)$, whichever places $p',q'$ in largest-to-smallest order. In this case, we have $p-2q=9$ and thus the ancestor pair is $(p',q')=(9,2)$ and therefore the ancestor triple is $(p'^2-q'^2,2p'q',p'^2+q'^2)=(77,36,85)$.
The only thing I don't [understand] is the last part where it's given as "One may thus equate numerators with numerators and denominators with denominators, giving Euclid's formula"
It uses unique fractionization $\Rightarrow$ uniqueness of reduced fractions (with denominators $> 0),\,$ i.e.
$\qquad\qquad \begin{align}\gcd(\color{#c00}{c,b})=1\\ \gcd(j,k)= 1\end{align}$, $\ \ \dfrac{c}b = \dfrac{j}k\ \Rightarrow\ \begin{align} c&\,=\,j\\ b &\,=\, k\end{align},\ \ \ {\rm for}\ \ b,c,j,k\in \Bbb Z,\ b,k > 0$
Follow the link for a simple proof using Euclid's Lemma (hint: $\,j = ck/b\,\Rightarrow\,\color{#c00}{b\mid c}\,k\,\Rightarrow\,b\mid k)$
Remark $ $ A more conceptual way to derive this parametrization of Pythagorean triples is to employ arithmetic of Gaussian integers $\,\Bbb Z[i] = \{ a + b\,i\,: a,b\in\Bbb Z\}$. Like integers they enjoy (Euclidean) division (with smaller remainder) and this implies they too satisfy the analog of the Fundamental Theorem of Arithmetic = existence and uniqueness of factorization into primes (= irreducibles). This implies that coprime factors of a square must themselves be squares (up to unit factors $\,\pm1,\pm i)$
Thus if $\ z^2 = x^2 + y^2 = (x-y\,i) (x+ y\,i) $ and $\,x,y\,$ are coprime then one easily checks that $\,x-y\,i,\,x+y\,i\,$ are coprime, so being coprime factors of the square $\,z^2$ they must themselves be squares (up to a unit factor). Thus e.g. $\ x + y\ i\, =\, (m + n\ i)^2 =\ m^2 - n^2 + 2mn\, i,\,$ hence $\,x = m^2-n^2,\ y = 2mn\,$ (using the unit factor $1$; using the other unit factors $\, -1,\pm i\,$ merely changes signs or swaps $\,x,y\,$ values). Notice how very simple the solution is from this perspective.
This is a simple prototypical (arithmetical) example of the simplification that results from the transformation of nonlinear problems into linear problems by working in larger algebraic extension rings. See here for some further discussion of such.
Best Answer
It doesn't seem cumbersome to me. Just loop on $m$ starting at $2$. Loop on $n$ starting at $1$ or $2$ depending on the parity of $m$ and going up to $m-1$. Check the GCD using Euclid's algorithm. If it is $1$ you have a primitive set, so calculate $a,b,c$. If you want all sets up to some $N$, multiply $a,b,c$ by all values up to $\lfloor \frac Nc \rfloor$. Stop the $m$ loop at $\sqrt N$