An infinite family of generalized Pell equations

Question

Let $m,n,x,y,z$ be positive integers such that $y=m(z-1)+1$ and $yz-1=n(y^2-x^2)$.

If you fix $m$ and $n$ then this is a generalized Pell equation in $x$ and $z$ (i.e., a quadratic Diophantine equation of the sort that can be solved by this website).

---

I have 2 specific conjectures, which I have checked computationally for $m,n\leq100$, $z\leq10^7$.

Conjecture 1: If $n\geq2$ then $x,y,z$ are all odd.

Conjecture 2: If $(m,n)\neq(1,1)$ then $z$ is indivisible by $4$.

A proof or counterexample of either of these conjectures would be much appreciated.

---

The table below gives the first four $(x,y,z)$ solutions for fixed $(m,n)$:

$$\begin{array}{c|c|c|c}&n=1&n=2&n=3\\\hline m=1& \begin{array}{c}(1,1,1)\\(1,2,2)\\(1,3,3)\\(1,4,4)\end{array}& \begin{array}{c}(1,1,1)\\(5,7,7)\\(29,41,41)\\(169,239,239)\end{array}& \begin{array}{c}(1,1,1)\\(9,11,11)\\(89,109,109)\\(881,1079,1079)\end{array}\\\hline m=2& \begin{array}{c}(1,1,1)\\(2,3,2)\\(23,33,17)\\(64,91,46)\end{array}& \begin{array}{c}(1,1,1)\\(167,193,97)\\(32397,37409,18705)\\(6284851,7257121,3628561)\end{array}& \begin{array}{c}(1,1,1)\\(439,481,241)\\(211597,231793,115897)\\(101989315,111723697,55861849)\end{array}\\\hline m=3& \begin{array}{c}(1,1,1)\\(3,4,2)\\(69,85,29)\\(287,352,118)\end{array}& \begin{array}{c}(1,1,1)\\(417,457,153)\\(200993,220177,73393)\\(96878209,106124761,35374921)\end{array}& \begin{array}{c}(1,1,1)\\(1053,1117,373)\\(1215161,1288873,429625)\\(1402294741,1487358181,495786061)\end{array} \end{array}$$

In the comments, Will Jagy observed that if $n\geq2$ then the $x$ values satisfy a linear recurrence $x_{n+2}=Cx_{n+1}-x_n$ (and you can do something similar for $y$ and $z$). The table below gives the values of $C$:

$$\begin{array}{c|c|c|c|c|c|c|c|c|c}&n=2&n=3&n=4&n=5&n=6&n=7&n=8&n=9&n=10\\\hline m=1&6&10&14&18&22&26&30&34&38\\\hline m=2&194&482&898&1442&2114&2914&3842&4898&6082\\\hline m=3&482&1154&2114&3362&4898&6722&8834&11234&13922\\\hline m=4&898&2114&3842&6082&8834&12098&15874&20162&24962\\\hline m=5&1442&3362&6082&9602&13922&19042&24962&31682&39202\\\hline m=6&2114&4898&8834&13922&20162&27554&36098&45794&\\\hline m=7&2914&6722&12098&19042&27554&37634&49282&&\\\hline m=8&3842&8834&15874&24962&36098&49282&&&\\\hline m=9&4898&11234&20162&31682&45794&&&&\\\hline m=10&6082&13922&24962&39202&&&&&\end{array}$$

Notice the symmetry!

Gottfried Helms observed the formula $C=4(2mn-1)^2-2$ for $m,n\geq2$.

When $m=1$, it looks like $C=4n-2$ for $n\geq2$.

---

This answer proves the $m=1$ case of conjecture 2, using methods from Diophantine approximation.

Answer 1

Phew! Thanks to the help of Will Jagy and Gottfried Helms, I have a complete solution!

Fix positive integers $m$ and $n$. We will study the equation \begin{equation} \tag{1} mn(y-x)(y+x)=(y-1)(y+m). \end{equation} Classification of Small Solutions

We begin by classifying the integer solutions to (1) with small $y$ coordinate.

Case 1: $y>1$

Assume that $y>1$. Then the right hand side of (1) is positive, so the left hand side of (1) is also positive. Then $y^2-x^2$ is positive, so $y^2-x^2\geq2y-1$. This gives $$(y-1)(y+m)=mn(y^2-x^2)\geq mn(2y-1).$$ Dividing through by $2y-1$ gives $$\frac{(y-1)(y+m)}{2y-1}\geq mn.$$ Let $f(y)=(y-1)(y+m)/(2y-1)$. Note that $f$ is increasing on $(\frac{1}{2},\infty)$. Let $y_0=2mn-m$. Then $$f(y_0)=\frac{2mn(y_0-1)}{2y_0-1}<mn\leq f(y)$$ so $y_0<y$. Thus, $y\geq2mn-m+1$.

Case 2: $-m<y<1$

Assume that $-m<y<1$. Then the right hand side of (1) is negative, so the left hand side of (1) is also negative. Then $y^2-x^2$ is negative, so $y^2-x^2\leq2y-1$. This gives $$(y-1)(y+m)=mn(y^2-x^2)\leq mn(2y-1).$$ Dividing through by $2y-1$ gives $$\frac{(y-1)(y+m)}{2y-1}\geq mn.$$ Let $f(y)=(y-1)(y+m)/(2y-1)$. Note that $f$ is increasing on $(-\infty,\frac{1}{2})$. Let $y_0=0$. Then $$f(y_0)=m\leq mn\leq f(y)$$ so $y_0\leq y$. Thus, $y=0$.

Case 3: $y<-m$

Assume that $y<-m$. Then the right hand side of (1) is positive, so the left hand side of (1) is also positive. Then $y^2-x^2$ is positive, so $y^2-x^2\geq-2y-1$. This gives $$(y-1)(y+m)=mn(y^2-x^2)\geq mn(-2y-1).$$ Dividing through by $-2y-1$ gives $$\frac{(y-1)(y+m)}{-2y-1}\geq mn.$$ Let $f(y)=(y-1)(y+m)/(-2y-1)$. Note that $f$ is decreasing on $(-\infty,-\frac{1}{2})$. Let $y_0=-2mn-m+2$. Then $$f(y_0)=\frac{(y_0-1)(y_0+m)}{-2y_0-1}\leq\frac{(y_0+1)(y_0+m-2)}{-2y_0-1}=\frac{-2mn(y_0+1)}{-2y_0-1}<mn\leq f(y)$$ so $y_0>y$. Thus, $y\leq-2mn-m+1$.

Completing the Classification

Theorem 1: If $(x,y)\in\mathbb{Z}^2$ satisfies (1) then exactly one of the following holds:

• $y\geq2mn-m+2$,
• $(x,y)=(\pm m,m+1)$ and $n=1$,
• $(x,y)=(\pm1,1)$,
• $(x,y)=(\pm1,0)$ and $n=1$,
• $(x,y)=(\pm m,-m)$,
• $(x,y)=(\pm(3m-2),-3m+1)$ and $n=1$,
• $y\leq-2mn-m$.

Proof: We have already shown that either $y\geq2mn-m+1$ or $y=1$ or $y=0$ or $y=-m$ or $y\leq-2mn-m+1$.

• Assume that $y=2mn-m+1$. Then (1) becomes $$mn(y-x)(y+x)=(2mn-m)(2mn+1).$$ Dividing through by $m$ and reducing modulo $n$ forces $n=1$. Then $y=m+1$ and solving for $x$ gives $x=\pm m$.
• If $y=1$ then solving for $x$ gives $x=\pm1$.
• If $y=0$ then solving for $x$ gives $nx^2=1$. This forces $n=1$ and $x=\pm1$.
• If $y=-m$ then solving for $x$ gives $x=\pm m$.
• Assume that $y=-2mn-m+1$. Then (1) becomes $$mn(y-x)(y+x)=(2mn+m)(2mn-1).$$ Dividing through by $m$ and reducing modulo $n$ forces $n=1$. Then $y=-3m+1$ and solving for $x$ gives $x=\pm(3m-2)$. $\blacksquare$

Mapping Solutions

We now define functions $f,g\colon\mathbb{R}^2\to\mathbb{R}^2$ that preserve integer solutions to (1). Our goal is to show that every integer solution to (1) can be obtained via $f$ and $g$ from one of the five solutions listed in Theorem 1.

Let $R=2mn-1$. Let $f,g\colon\mathbb{R}^2\to\mathbb{R}^2$ be defined by $$f\left(\begin{bmatrix}x\\y\end{bmatrix}\right)=\begin{bmatrix}(2R^2-1)x+2R(R-1)y-2(m-1)R\\(2R^2-1)y+2R(R+1)x-2(m-1)(R+1)\end{bmatrix},$$ $$g\left(\begin{bmatrix}x\\y\end{bmatrix}\right)=\begin{bmatrix}(2R^2-1)x-2R(R-1)y+2(m-1)R\\(2R^2-1)y-2R(R+1)x-2(m-1)(R+1)\end{bmatrix}.$$

Lemma 2: The functions $f$ and $g$ preserve solutions to (1).

Proof: The algebraic manipulations are rather involved, so I instead wrote up a proof in Lean 3.

import tactic.ring

variables {A : Type*} [comm_ring A]

lemma lem0 {a b c d : A} (h : a = b) : c - a = d - b → c = d :=
by { rw h, exact sub_left_inj.mp }

structure sol (m n : A) :=
(x y : A)
(h : m * n * (y - x) * (y + x) = (y - 1) * (y + m))

def R (m n : A) : A := 2 * m * n - 1

def sol.forward {m n : A} (s : sol m n) : sol m n :=
{ x := (2*(R m n)*(R m n)-1)*s.x+2*(R m n)*((R m n)-1)*s.y-2*(m-1)*(R m n),
  y := (2*(R m n)*(R m n)-1)*s.y+2*(R m n)*((R m n)+1)*s.x-2*(m-1)*((R m n)+1),
  h := by { dsimp only [R], apply lem0 s.h, ring } }

def sol.backward {m n : A} (s : sol m n) : sol m n :=
{ x := (2*(R m n)*(R m n)-1)*s.x-2*(R m n)*((R m n)-1)*s.y+2*(m-1)*(R m n),
  y := (2*(R m n)*(R m n)-1)*s.y-2*(R m n)*((R m n)+1)*s.x-2*(m-1)*((R m n)+1),
  h := by { dsimp only [R], apply lem0 s.h, ring } }

You just need to check that I correctly copied (1) (sol.h), $f$ (sol.forward), and $g$ (sol.backward). $\blacksquare$

We now give a few lemmas describing the relationship between $f$ and $g$.

Lemma 3: Let $\overline{(x,y)}=(-x,y)$. Then $f(\overline{p})=\overline{g(p)}$ for all points $p\in\mathbb{R}^2$.

Proof: This is clear from the definitions of $f$ and $g$. $\blacksquare$

Lemma 4: The functions $f$ and $g$ are inverses to each other.

Proof: We can write $f$ in terms of matrix multiplication as $$f\left(\begin{bmatrix}x\\y\end{bmatrix}\right)=\begin{bmatrix}2R^2-1&2R(R-1)\\2R(R+1)&2R^2-1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}-2(m-1)\begin{bmatrix}R\\R+1\end{bmatrix}.$$ Note that the $2\times2$ matrix has determinant $(2R^2-1)^2-4R^2(R^2-1)=1$. Then the inverse function to $f$ is given by \begin{align*} f^{-1}\left(\begin{bmatrix}x\\y\end{bmatrix}\right)&=\begin{bmatrix}2R^2-1&2R(R-1)\\2R(R+1)&2R^2-1\end{bmatrix}^{-1}\left(\begin{bmatrix}x\\y\end{bmatrix}+2(m-1)\begin{bmatrix}R\\R+1\end{bmatrix}\right)\\ &=\begin{bmatrix}2R^2-1&-2R(R-1)\\-2R(R+1)&2R^2-1\end{bmatrix}\left(\begin{bmatrix}x\\y\end{bmatrix}+2(m-1)\begin{bmatrix}R\\R+1\end{bmatrix}\right)\\ &=\begin{bmatrix}2R^2-1&-2R(R-1)\\-2R(R+1)&2R^2-1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}+2(m-1)\begin{bmatrix}2R^2-1&-2R(R-1)\\-2R(R+1)&2R^2-1\end{bmatrix}\begin{bmatrix}R\\R+1\end{bmatrix}\\ &=\begin{bmatrix}2R^2-1&-2R(R-1)\\-2R(R+1)&2R^2-1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}+2(m-1)\begin{bmatrix}R\\-(R+1)\end{bmatrix}=g\left(\begin{bmatrix}x\\y\end{bmatrix}\right). \end{align*} Thus, $f$ and $g$ are inverses to each other. $\blacksquare$

By Lemma 4, fixed points of $f$ are fixed points of $g$, and vice versa.

Lemma 5: If $(x,y)$ is a fixed point of $f$ and $g$ then $x=0$.

Proof: Suppose that $(x,y)$ is a fixed point of $f$ and $g$. Then $$\begin{bmatrix}2R^2-1&2R(R-1)\\2R(R+1)&2R^2-1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}-2(m-1)\begin{bmatrix}R\\R+1\end{bmatrix}=f\left(\begin{bmatrix}x\\y\end{bmatrix}\right)=\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}1&0\\0&1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}.$$ Rearranging terms and dividing through by 2 gives \begin{equation}\tag{2}\begin{bmatrix}R^2-1&R(R-1)\\R(R+1)&R^2-1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}=(m-1)\begin{bmatrix}R\\R+1\end{bmatrix}.\end{equation} Note that the $2\times2$ matrix has determinant $(R^2-1)^2-R^2(R^2-1)=1-R^2$. If $R=1$ then $(m,n)=(1,1)$, so (2) becomes $$\begin{bmatrix}0&0\\2&0\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}0\\0\end{bmatrix}.$$ This forces $x=0$. Now suppose that $R>1$. Then the $2\times2$ matrix in (2) is invertible, so (2) has a unique solution. In other words, $f$ and $g$ have a unique fixed point $p$. However, Lemma 3 shows that $\overline{p}$ is also a fixed point of $f$ and $g$. Then $p=\overline{p}$, so $x=0$. Alternatively, you can solve (2) explicitly and get $(x,y)=(0,(m-1)/(R-1))$. $\blacksquare$

Lemma 6: The functions $f,g\colon\mathbb{R}^2\to\mathbb{R}^2$ have no fixed points satisfying equation (1).

Proof: Let $(x,y)$ be a fixed point of $f$ and $g$. By Lemma 5, $x=0$. However, setting $x=0$ in (1) gives $$mny^2=(y-1)(y+m).$$ The quadratic $(mn-1)y^2-(m-1)y+m$ has discriminant $(m-1)^2-4m(mn-1)<0$. $\blacksquare$

We can now prove the main theorem.

Theorem 7: Any integer solution to (1) can be mapped via $f$ and $g$ to one of the five solutions listed in Theorem 1.

Proof: Since the quadratic $(mn-1)y^2-(m-1)y+m$ has negative discriminant, the graph of (1) consists of the two curves $x=\pm\sqrt{y^2-((y-1)(y+m)/(mn))}$. Since $f(-1,1)=(1-2mR,1-2m(R+1))$ has negative $x$-coordinate, $f$ and $g$ do not swap the two curves. By Lemma 6, $f$ and $g$ do not flip either of the two curves, and each shift each curve either up or down. Since $f(-1,1)=(1-2mR,1-2m(R+1))$ has negative $y$-coordinate, $f$ shifts the left curve downward. Then by Lemma 4, $g$ shifts the left curve upward. Then by Lemma 3, $f$ shifts the right curve upward and $g$ shifts the right curve downward.

It remains to show that $f$ and $g$ can't jump over the gap $-2mn-m+1\leq y\leq2mn-m+1$. By Lemmas 3 and 4, it suffices to show that if $y_0=2mn-m+1$ and $x_0=\sqrt{y_0^2-((y_0-1)(y_0+m)/(mn))}$ then the $y$-coordinate of $g(x_0,y_0)$ is at least $-2mn-m+1$. In other words, we need to prove the inequality $$(2R^2-1)y_0-2R(R+1)x_0-2(m-1)(R+1)\geq-2mn-m+1.$$ Rearranging terms gives the equivalent inequality $$x_0\leq\frac{(2R^2-1)y_0-2(m-1)(R+1)+2mn+m-1}{2R(R+1)}.$$ Simplifying the right hand side gives the equivalent inequality $$x_0\leq2mn-m.$$ To prove this, we directly compute $$y_0^2-x_0^2=\frac{(y_0-1)(y_0+m)}{mn}=\frac{(2mn-m)(2mn+1)}{mn}=4mn-2m+2-\frac{1}{n}\geq4mn-2m+1=y_0^2-(y_0-1)^2.$$ Then $x_0\leq y_0-1$ as desired. $\blacksquare$

Consequences

We can combine Theorems 1 and 7 to make conclusions about all integer solutions to (1).

Theorem 8: Let $(x,y)\in\mathbb{Z}^2$ satisfy (1).

• $y\equiv0,1\pmod{m}$.
• If $n\neq1$ and $y\equiv0\pmod{m}$ then $x\equiv y\equiv m\pmod{2}$.
• If $n\neq1$ and $y\equiv1\pmod{m}$ then $x\equiv y\equiv1\pmod{2}$.
• Either way, if $n\neq1$ then $x\equiv y\pmod{2}$.

Proof: Theorem 7 states that any integer solution to (1) can be mapped via $f$ and $g$ to one of the five solutions listed in Theorem 1. Note that $f$ and $g$ preserve the value of $y\pmod{m}$, $x\pmod{2}$, and $y\pmod{2}$. Since each of the five solutions listed in Theorem 1 satisfy the theorem statement, all integer solutions to (1) satisfy the theorem statement. $\blacksquare$

Theorem 9: Let $m,n,x,y,z$ be positive integers such that $y=m(z-1)+1$ and $yz-1=n(y^2-x^2)$.

• If $n\geq2$ then $x,y,z$ are all odd.
• If $(m,n)\neq(1,1)$ then $z$ is indivisible by 4.

Proof: Note that $$mn(y^2-x^2)=m(yz-1)=y(mz)-m=y(y+m-1)-m=(y-1)(y+m)$$ so our previous work applies. Furthermore, we are in the $y\equiv1\pmod{m}$ case. By Theorem 8, if $n\geq2$ then $x\equiv y\equiv1\pmod{2}$. Then the equation $yz-1=n(y^2-x^2)$ forces $z\equiv1\pmod{2}$.

It remains to show that if $m\geq2$ and $n=1$ then $z$ is indivisible by 4. We know that $(x,y)$ can be mapped via $f$ and $g$ to one of the five solutions listed in Theorem 1. Since we are in the $y\equiv1\pmod{m}$ case, and since $m\geq2$, we know that $(x,y)$ can be mapped via $f$ and $g$ to $(\pm m,m+1)$ or $(\pm1,1)$ or $(\pm(3m-2),-3m+1)$. Since $(x,y)$ is in the first quadrant, we know that $(x,y)$ can be obtained from $(1,1)$ or $(m,m+1)$ via repeated application of $f$ (we can discard $(3m-2,-3m+1)$ since $g(m,m+1)=(3m-2,-3m+1)$). The solution $(1,1)$ has $z=1$. The solution $(m,m+1)$ has $z=2$.

We now check how applying $f$ changes $z$. Applying $f$ has the effect $$y\mapsto(2R^2-1)y+2R(R+1)x-2(m-1)(R+1).$$ In terms of $z$, this is $$z\mapsto\frac{(2R^2-1)(m(z-1)+1)+2R(R+1)x-2(m-1)(R+1)-1}{m}+1$$ which simplifies to $$z\mapsto(2R^2-1)z+4nRx-4(m-1)nR.$$ In particular, the value of $z\pmod{4}$ does not change. $\blacksquare$

Answer 2

This is not an answer, just a comment to supply more data.

I used Pari/GP to detect the polynomial-structure in your C-table, and find the following table containing additional info for the column $n=1$:

  m \ n|    1     2     3      4      5      6      7      8  |
    -  +    -     -     -      -      -      -      -      -  +
    1  |    2    34    98    194    322    482    674    898  |   // ????? see note
    2  |   34   194   482    898   1442   2114   2914   3842  |
    3  |   98   482  1154   2114   3362   4898   6722   8834  |
    4  |  194   898  2114   3842   6082   8834  12098  15874  |
    5  |  322  1442  3362   6082   9602  13922  19042  24962  |
    6  |  482  2114  4898   8834  13922  20162  27554  36098  |
    7  |  674  2914  6722  12098  19042  27554  37634  49282  |
    8  |  898  3842  8834  15874  24962  36098  49282  64514  |
    -  +  ...   ...   ...    ...    ...    ...    ...    ...  +
          ???

This table is produced by the small Pari/GP function

Cpol(m,n='x)= (4*m*n)^2-16*m*n+2 \\ reproduces OP's table for m,n>=2
or, a bit more symmetric,
Cpol(m,n='x)= 4*(2*m*n-1)^2 - 2 \\ reproduces OP's table for m,n>=2
$\qquad \qquad$ _{(leaving n in the function-call indeterminate gives the $m$'th polynomial in $x$)}

Surprisingly, the first row, (for $m=1$) does not agree with your first row; and the same might occur, if you would supply a first column with $n=1$ of your data (however looking numerically for solutions when $n=1$ agrees with the first column here).
I'm checking further to see what's going on with this.

---

update1
I have looked at the solutions for $n=1$ and $m \gt 1$ .

First, with a chosen $m$, I take, from sequential search, the first 6 solutions for $(x,y,z)$. We find, that assuming we have 2 interlaced subsequences, that subsequences $x_a$ and $x_b$ have a simple recursion rule; for $x_a$ and $x_b$ we have $$x_{a,k+2}=C_{m,1}\cdot x_{a,k+1} - 1 \cdot x_{a,k} \\ x_{b,k+2}=C_{m,1}\cdot x_{b,k+1} - 1 \cdot x_{b,k} \\\tag {1.x}$$

The empirical data agree with the $C_{m,n}$-values from my proposal.

For the interlaced $y$- and $z$- sequences we have to include a constant value $D_{y,m}$ and $D_{z,m}$ in the recursion formula: $$ y_{a,k+2}=C_{m,1}\cdot y_{a,k+1} - 1 \cdot y_{a,k} + D_{y,m}\\ y_{b,k+2}=C_{m,1}\cdot y_{b,k+1} - 1 \cdot y_{b,k} + D_{y,m}\\ \tag {1.y}$$ and $$ z_{a,k+2}=C_{m,1}\cdot z_{a,k+1} - 1 \cdot z_{a,k} + D_{z,m}\\ z_{b,k+2}=C_{m,1}\cdot z_{b,k+1} - 1 \cdot z_{b,k} + D_{z,m}\\ \tag {1.z}$$

The complete recursion can be done in a matrix-formula, where the recursion-depth can be formalized by an initial matrix $I$ and the multiplicatormatrix $M$ such that $$ S_{h+1} = I \cdot M^h $$ gives the $h$'th iteration.

Solutions:

• Let $m=2$. The first few empirically found solutions for $(x,y,z)$ are empirically the following sequences, obviously by 2 recursive intertwined subsequences:

  x=  1  2  23  64   781  2174  26531   73852
  y=  1  3  33  91  1105  3075  37521  104443
  z=  1  2  17  46   553  1538  18761   52222
  

The matrix $M$ is $$ M = \small{ \begin{bmatrix} 1&0&1 \\ 0&0&-1\\0&1&C_{m,1}\end{bmatrix} = \begin{bmatrix} 1&0&1 \\ 0&0&-1\\0&1&34\end{bmatrix} }$$

The initmatrix $I$ is composed of the first two empirical values $[x_a,y_a,z_a]$ and $[x_b,y_b,z_b]$, and additionally prepended by the -at the moment- unknown constants $D$ taken as zero: $$ I_{\text{init}}=\small \begin{bmatrix} 0&x_{a,1}&x_{a,2} \\ D_y&y_{a,1}&y_{a,2}\\D_z&z_{a,1}&z_{a,2} \\ \hline 0&x_{b,1}&x_{b,2} \\ D_y&y_{b,1}&y_{b,2}\\D_z&z_{b,1}&z_{b,2} \end{bmatrix} = \begin{bmatrix} 0&1&23\\0&1&33\\0&1&17 \\ \hline 0&2&64\\ 0&3&91\\0&2&46 \end{bmatrix} $$ If we compute $S_2 = I \cdot M$ we get a wrong result for $y_{a,3},z_{a,3},y_{b,3},z_{a,3}$ - but this differences are then simply the new constant values for $D$ and this gives then $$ \small I= \begin{bmatrix} 0&1&23\\-16&1&33\\-24&1&17 \\ \hline 0&2&64\\ -16&3&91\\-24&2&46 \end{bmatrix}$$

With this we can compute the set of solutions using $S_{h+1}=I \cdot M^h$ for $h=0..7$ (or arbitrarily more); I show this having the two subsets as submatrices: $$\small \begin{bmatrix} o: & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & \cdots \\ \hline x_a: & 1 & 23 & 781 & 26531 & 901273 & 30616751 & 1040068261 & 35331704123 & \cdots \\ y_a: & 1 & 33 & 1105 & 37521 & 1274593 & 43298625 & 1470878641 & 49966575153 & \cdots \\ z_a: & 1 & 17 & 553 & 18761 & 637297 & 21649313 & 735439321 & 24983287577 & \cdots \\ x_b: & 2 & 64 & 2174 & 73852 & 2508794 & 85225144 & 2895146102 & 98349742324 & \cdots \\ y_b: & 3 & 91 & 3075 & 104443 & 3547971 & 120526555 & 4094354883 & 139087539451 & \cdots \\ z_b: & 2 & 46 & 1538 & 52222 & 1773986 & 60263278 & 2047177442 & 69543769726 & \cdots \end{bmatrix}$$

Note that we could even insert negative values for $h$ to continue the matrix of solutions to the left side.

---

• Let $m=3$. The first few empirically found solutions for $(x,y,z)$ are empirically the following sequences, obviously again by 2 recursive intertwined subsequences:

  x=    1  3  69  287  6761  28123
  y=    1  4  85  352  8281  34444
  z=    1  2  29  118  2761  11482
  

The matrix $M$ is $$ M = \begin{bmatrix} 1&0&1 \\ 0&0&-1\\0&1&98\end{bmatrix} $$

The initmatrix $I$ is made as before and this gives then $$ I = \small \begin{bmatrix} 0 & 1 & 69 \\ -48 & 1 & 85 \\ -80 & 1 & 29 \\ \hline 0 & 3 & 287 \\ -48 & 4 & 352 \\ -80 & 2 & 118 \end{bmatrix}$$

With this we can compute the set of solutions using $S_{h+1}=I \cdot M^h$ for $h=0..7$ (or arbitrarily more); again I show this having the two subsets as submatrices: $$\small \begin{bmatrix} o: & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & \cdots \\ \hline x_a: & 1 & 69 & 6761 & 662509 & 64919121 & 6361411349 & 623353393081 & 61082271110589 \\ y_a: & 1 & 85 & 8281 & 811405 & 79509361 & 7791105925 & 763448871241 & 74810198275645 \\ z_a: & 1 & 29 & 2761 & 270469 & 26503121 & 2597035309 & 254482957081 & 24936732758549 \\ \hline x_b: & 3 & 287 & 28123 & 2755767 & 270037043 & 26460874447 & 2592895658763 & 254077313684327 \\ y_b: & 4 & 352 & 34444 & 3375112 & 330726484 & 32407820272 & 3175635660124 & 311179886871832 \\ z_b: &2 & 118 & 11482 & 1125038 & 110242162 & 10802606758 & 1058545220042 & 103726628957278 \end{bmatrix} $$