Polynomial parametrization of the solutions to $yz=x^2+x\pm 1$ If a Diophantine equation has infinitely many integer solutions, how to describe them all? One standard approach is polynomial parametrization. For example, all integer solutions to the equation
$$
yz=x^2
$$
are given by $x=uvw$, $y=uv^2$, $z=uw^2$ for some integers $u,v,w$. More formally, a subset $S \subset {\mathbb Z}^n$ is a polynomial family if there exists polynomials $P_1,\dots,P_n$ in some variables $u_1,\dots,u_k$ and integer coefficients such that $(x_1,\dots,x_n) \in S$ if and only if there exists integers $u_1,\dots,u_k$ such that $x_i=P_i(u_1,\dots,u_k)$ for $i=1,\dots,n$.
For some equations the solution set is not a polynomial family but is a finite union of polynomial families. The simplest example is the equation $xy=0$ with solutions $(x,y)=(u,0)$ and $(0,u)$.
In 2010, Vaserstein https://annals.math.princeton.edu/wp-content/uploads/annals-v171-n2-p07-s.pdf answered a long-standing open question and showed that the solution set to the equation
$$
xy - zt = 1
$$
is a polynomial family. As a corollary, he showed that the solutions to many equations, including $yz=x^2+a$ for any $a$ and $xy-zt=a$ for any $a$, are the finite unions of polynomial families. As noted by Fedor Pertov in the comment, this result also implies parametrization of the solution set of the equation
$$
yz=x^2+x.
$$
The simplest examples that seems to be not explicitly covered by this paper are equations
$$
yz=x^2+x+1
$$
and
$$
yz=x^2+x-1.
$$
The question is, for each of these equations, whether its solution set is a finite union of polynomial families? Or, in simple words, can we write down all solutions using polynomial expressions with parameters?
 A: For $x^2+x+1=yz$ we may factorize LHS in the unique factorization domain $\mathbb{Z}[\omega]$, where $\omega=e^{2\pi i/3}$: $$(x-\omega)(x-\omega^2)=yz.$$
Denote by $A$ the greatest common divisor of $x-\omega$ and $y$, say $x-\omega=AB$, $y=AC$. Then $B$ and $C$ are coprime, and $B(x-\omega^2)=Cz$. Thus $C$ divides $x-\omega^2$, that is, $x-\omega^2=CD$, $z=BD$. Since $AC$ is real, we get $A/\overline{C}\in \mathbb{R}$. Since both $A$ and $\overline{C}$ belong to $\mathbb{Z}[\omega]$, the line through 0 and $A$ intersects $\mathbb{Z}[\omega]$ by an infinite cyclic subgroup. So we get $A=pT$, $\overline{C}=qT$ (equivalently, $C=q\overline{T}$) for certain $p,q\in \mathbb{Z}$ and $T\in \mathbb{Z}[\omega]$. Thus $p$ divides $x-\omega$, but $(x-\omega)/p\in \mathbb{Z}[\omega]$ only if $p=\pm 1$, analogously $q=\pm 1$. So, $C=\pm \overline{A}$, analogously $D=\pm \overline{B}$. So, we should parametrize the solutions of $$x-\omega=AB,$$
the rest is automatic. If $A=u+v\omega$, $B=u_1+v_1\omega^2$, then $$AB=(u+v\omega)(u_1+v_1\omega^2)=uu_1+vv_1+(vu_1-uv_1)\omega-uv_1.$$
It equals $x-\omega$ if and only if $vu_1-uv_1=-1$ (and $x=uu_1+vv_1-uv_1$), which may be parametrized by the cited theorem of Vaserstein..
A: I will prove that, more generally, for any integer parameters $a,b,c,d$, the solution set to any equation in the form
$$
ax^2 + bx + c = dyz
$$
is a finite union of polynomial families.
Lemma: Let $S \subset Z^n$ be a finite union of polynomial families. Fix integers $a_1,\dots,a_n$, and positive integers $b_1,\dots,b_n$. Let $S' \subset S$ be the set of all $x=(x_1,\dots, x_n) \in S$ such that $x_i \equiv a_i (\text{mod }b_i)$ for all $i$. Then $S'$ is also a finite union of polynomial families.
Proof: We may assume that $S$ is a polynomial family. By definition, there exists polynomials $P_1,\dots,P_n$ in some variables $u=(u_1,\dots,u_k)$ such that, for every $x \in S$, we have $x_i=P_i(u), i=1,\dots,k$ for some $u \in {\mathbb Z}^k$. Let $m$ be the least common multiple of $b_1,\dots,b_n$. Let $Q$ be the set of vectors $q=(q_1,\dots q_k)$ such that all $q_i$ are integers such that $0 \leq q_i < m$. Let $Q' \subset Q$ be the set of $q \in Q$ such that $P_i(q_1,\dots,q_k)\equiv a_i (\text{mod } b_i)$ for every $i$. For each $q \in Q'$, do a substitution $u_i \to mv_i + q_i$ where $v_i$ is a new variable. Let $P^q_i(v_1,\dots,v_k)=P_i(mv_1 + q_1, \dots, mv_k+q_k)$. This way we have represented $S'$ as the union of $|Q'|$ polynomial families defined by $(P^q_1, \dots, P^q_n)$.
Corollary If the solution set $S$ to the equation $P(x_1,\dots,x_n)$ is a finite union of polynomial families, then so is the solution set $S'$ to the equation $P(b_1x_1-a_1,\dots,b_nx_n-a_n)=0$ for integers $a_1,\dots,a_n$ and non-zero integers $b_1,\dots,b_n$.
Now multiply the equation $ax^2+bx+c = dyz$ by $4a$ to get $4adyz = (2ax+b)^2 + 4ac-b^2$. This is the equation of the form $yz=x^2-D$ with $D=b^2-4ac$ with the restriction that $x \equiv -b (\text{mod } 2a)$ and $y \equiv 0 (\text{mod }4ad)$. As noted in the question, Vaserstein proved that the solution set to $yz=x^2-D$ is a finite union of polynomial families. It is left to apply the corollary.
