Looking through some of my posts recently I found myself thinking more on certain research, which is necessary to show how easy it is to keep going to infinity if you have solutions to x2 - Dy2 = 1, where x = 1 mod D or x = -1 mod D.
So I'm echoing things already discussed on the blog, but pulling out particular things to highlight.
Given x2 - Dy2 = 1, and x = 1 mod D or x = -1 mod D:
(D-1)j2 + (j+1)2 = (x+y)2, where j = ((x+Dy)-1)/D
(D-1)j2 + (j-1)2 = (x+y)2, where j = ((x+Dy)+1)/D
which gives a Pythagorean triple when D-1 is a square. The condition that x = 1 mod D or x = -1 mod D is always true if D is prime or twice a prime. For other cases though, when that requirement is not met, the following analysis does not apply.
Starting with the first relation I note that:
(j+1)2 - (x+y)2 = 0 mod D-1
Now consider any prime factor p of D-1, of course it also follows then that:
(j+1)2 - (x+y)2 = 0 mod p
But now consider j = ((x+Dy)-1)/D, and since D-1 = 0 mod p, I have j = x+y-1 mod p, which forces all prime factors of D-1 into: (j+1 -(x+y)).
Going back to my original full expression:
(D-1)j2 + (j+1)2 = (x+y)2, which is, (D-1)j2 = (x+y - (j+1))(x+y + (j+1))
From the original equation itself I have x2 - y2 = 1 mod p, so x+y must be coprime to D-1. And knowing that D-1 is a factor of (x+y - (j+1)), it must be coprime to (x+y + j+1) for all its primes except 2. Also then (x+y + j+1) is forced to be a perfect square or twice a perfect square.
So I have (x+y + j+1) = n2 or 2n2, for some integer n.
And (x+y - (j+1)) = (D-1)m2 or 2(D-1)m2, for some integer m.
That applies if x = 1 mod D.
Going to the other possibility, which is x = -1 mod D:
x+y + 2nm + 1 = 2n2, and x+y - (2nm+1) = 2(D-1)m2
and again subtracting the last from the first gives: 2(2nm+1) = 2n2 - 2(D-1)m2.
So n2 - 2nm - 1 - (D-1)m2 = 0, which gives, (n-m)2 - Dm2 = 1, which is the original equation back again, and I can solve for n to get:
n = m +/- sqrt(Dm2 + 1)
Those are just 2 of 4 possibles though, as (D-1)j2 + (j-1)2 = (x+y)2, then for the x = -1 mod D case gives:
(x+y + j-1) = n2 or 2n2, for some integer n.
And (x+y - (j-1)) = (D-1)m2 or 2(D-1)m2, for some integer m.
Originally I went down this path thinking it might lead to an easy general way to solve the original equation, but didn't notice one which is kind of how it fell to the side.
But more recently I've noted that you can use solutions for x and y to get more solutions to another Diophantine equation.
For instance with D=61:
x = 1766319049, y = 226153980, so x = 1766319049 = -1 mod 61.
Then j = 255110030.
And taking the positives with some square roots though I'm sure you can do something with the negative values:
sqrt((x+y+j-1)/2) = 33523 = n
sqrt((x+y-(j-1))/(2(D-1))) = 3805 = m
Notice that j = 2nm.
And (n-m)2 - 61m2 = -1.
297182 - 61*38052 = -1
The above complements another post which explains more, where I mention that you can get infinite cycles but didn't show how you calculate n and m. So the point of this post is to show how they are calculated.
And you know it is kind of wild to have that from my own original research.
Of course the math doesn't care. The results have always been true. It was just a matter of someone finding them.