Context: I am coding a solver (in Java) for generalized Pell equations of the form $\displaystyle x^2 - Dy^2 = N$ using a paper of

John Robertson as a reference. The Lagrange-Matthews-Mollin algorithm requires a complete solution set to the congruence $\displaystyle z^2 \equiv D\ (mod\ |m|)$ for various m.

So, I've handled factoring |m| into prime powers $\displaystyle p^k$ and am OK with Tonelli Shanks and Hensel lifting when $\displaystyle D \not \equiv 0\ (mod\ p)$, but I'm stuck on the supposedly trivial case $\displaystyle D \equiv 0\ (mod\ p)$. Say for example we want to find all solutions to $\displaystyle x^2 \equiv 98\ (mod\ 7^3)$. It is easy to check that the solutions are x $\displaystyle \equiv$ 21, 28, 70, 77, 119, 126, 168, 175, 217, 224, 266, 273, 315, 322 (mod $\displaystyle 7^3$). I can see that the solution set is simply $\displaystyle x \equiv \pm 21\ (mod\ 49)$. Likewise, I have determined that the solution set to $\displaystyle x^2 \equiv 98\ (mod\ 7^4)$ is $\displaystyle x \equiv \pm 70\ (mod\ 7^3)$. What would be an efficient algorithm for the general case? And will prime p = 2 need to be treated separately from primes p > 2?

I have limited number theory background, so I'm probably missing something basic. The literature I've consulted seems to pass off such problems as trivial or easy or irrelevant, leaving them as exercises to the reader, either explicitly or by omission.

After I've completed this step, straightforward application of the Chinese Remainder Theorem will allow me to complete the overall problem.

Any help would be appreciated.