Solving for modular solutions to x2 - Dy2 = 1, is surprisingly easy, with an alternate form. Consider
(z-y)2 - Dy2 = 1
as then z2 - 2zy + y2 - Dy2 = 1, so I can group to find D-1 makes an appearance:
z2 - 2zy - (D-1)y2 = 1, so z2 - 2zy - 1 = (D-1)y2
And then solve for y modulo D-1, to find:
2y = z-1(z2 - 1) mod D-1
So what is z? Well it is a residue of D-1, and intriguingly, the implication is that for every such residue except the trivial z = 1 mod D-1, if that residue is allowed, since the modular inverse blocks cases where z and D-1 share prime factors, then there is a solution for the original equation, which mathematicians commonly call Pell's Equation. I like to call it the two conics equation, as it gives hyperbolas or ellipses depending on the sign of D, though for integer solutions only, D is positive so you're considering hyperbolas.
Surprisingly if you do a search on "Pell's Equation" it does not appear that mathematicians noticed that you can so easily solve for y modularly, but it's not always clear that something is previously unknown for that reason. And it turns out that knowing y mod D-1 does not mean you have an explicit solution, as it can be greater than D-1.
However, I have also shown this result pushes rules upon quadratic residue pairs, which allows you to count them. And that's just one thing that I noticed quickly. Who knows how many ways this result impacts number theory at this point, but that is enough to suggest it was previously unknown.
And since I gave the rules for fundamental solutions to x2 - Dy2 = 1, it seems likely to me that not having a way to prove those--though they may have guessed them--mathematicians did not know how to solve for y modularly either.
But it gets better as generalizing may pull in explicit solutions of any size.
Consider now adding yet another variable k.
(z-ky)2 - Dy2 = 1
Then we solve to find:
2y = (kz)-1(z2 - 1) mod D-k2
So you can just find a series of values for y modulo a succession of values by shifting k, and since D can be less than k2, you have as many as you wish. And then you just use the Chinese remainder theorem to get an ever larger value for y, but how do you pick?
Well I'm not sure as I just had this idea but the math may not care, and may give you a solution regardless of which residues you use, which just may not be the fundamental solution. But that is just a conjecture at this point.
One way of considering it is there are an infinity of y's that will work, so your arbitrary pick of residues may be like dialing up one of them.
But if true it is a way to solve the original equation trivially at will, and you can work back from a larger solution to the fundamental solution if you really want it.
Intriguingly, these results also push against quadratic residues here as well, and who knows what the full impact is.