Just recently I had a need to solve the following problem:
Given $$$a$$$ and $$$m$$$, find all $$$x$$$ such that: $$$x^2 \equiv a \pmod{m}$$$.
After some time searching the Internet and carefully analyzing each case, I got a solution that works for any $$$a$$$ and $$$m > 0$$$ in $$$O(\sqrt{m})$$$ time, or in $$$O(\log^2(m))$$$ time if $$$a$$$ and $$$m$$$ are coprime (if there is a factorization of $$$m$$$).
Below are the cases I considered and links to articles (or discussions) for algorithms:
- $$$x^2 \equiv a \pmod{p}$$$, $$$a \ne 0$$$, $$$p$$$ — odd prime (link).
- $$$x^2 \equiv a \pmod{p^k}$$$, $$$a \ne 0$$$, $$$k > 0$$$, $$$\gcd(a, p) = 1$$$, $$$p$$$ — odd prime (link).
- $$$x^2 \equiv a \pmod{2^k}$$$, $$$a \ne 0$$$, $$$k > 0$$$, $$$\gcd(a, 2) = 1$$$ (link).
- $$$x^2 \equiv 0 \pmod{p^k}$$$, $$$k > 0$$$, $$$p$$$ — prime. Note that in this case $$$x^2$$$ must be divisible by $$$p^k$$$. Then $$$x$$$ must be divisible by $$$p^{\lceil\frac{k}{2}\rceil}$$$. From this we obtain that $$$x = i \cdot p^{\lceil\frac{k}{2}\rceil}$$$, $$$0 \le i < p^{\lfloor\frac{k}{2}\rfloor}$$$ .
- $$$x^2 \equiv p^s \pmod{p^k}$$$, $$$k > 0$$$, $$$0 < s < k$$$, $$$p$$$ — prime. Solutions exist only for $$$s = 2 \cdot t$$$. Proved by contradiction. Then it is obvious that $$$p^t$$$ is one of the solutions ($$$(p^t)^2 = p^s$$$). Let's find the remaining solutions. Let $$$x = p^t + y$$$, then $$$x^2 = p^s + 2 \cdot p^t \cdot y + y^2$$$. Note that $$$2 \cdot p^t \cdot y + y^2$$$ must be divisible by $$$p^k$$$. In this case, $$$x = p^t + i \cdot p^{\max(\lceil\frac{k}{2}\rceil, k - t - (p = 2))}$$$.
- $$$x^2 \equiv q \cdot p^s \pmod{p^k}$$$, $$$k > 0$$$, $$$0 < s < k$$$, $$$\gcd(q, p) = 1$$$, $$$p$$$ — prime. Let $$$x_1$$$ be the solution for $$$x^2 \equiv q \pmod{p^k}$$$, and let $$$x_2$$$ be the solution for $$$x^2 \equiv p^s \pmod{p^k}$$$. Then $$$x = x_1 \cdot x_2$$$ will be a solution to the original problem.
The last three cases had to be deduced independently, so there are no links. I was too lazy to document the proofs in this post, although in general they are not very complicated.
The general solution is found by using the Chinese Remainder Theorem and reducing to one of the above cases.
I will also give a complete implementation in Python (a function that solves the general case — discrete_sqrt):
