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?



Polynomial parametrization of the solutions to $yz=x^2+x\pm 1$
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..


