Taking square roots modulo $2^N$ I was trying to solve $y^2 - y \equiv 16 \pmod{512}$ by completing the square.
Here is my solution.
\begin{align}
   y^2 - y &\equiv 16 \pmod{512} \\
   4y^2 - 4y + 1 &\equiv 65 \pmod{512} \\
  (2y-1)^2 &\equiv 65 \pmod{512} \\
   2y - 1 &\equiv \pm 33 \pmod{512} &\text{Found by pointwise search.}\\
   2y &\equiv 34, -32 \pmod{512} \\
   y &\equiv 17, -16 \pmod{256} \\
   y &\in \{17, 273, 240, 496\} \pmod{512} &\text{These values need to be verified.}\\
   y &\in \{240, 273\}\pmod{512}
\end{align}
I had to solve $x^2 \equiv 65 \pmod{2^9}$ by a pointwise search. Is there any systematic method for solving equivalences of the form $x^2 \equiv a \pmod{2^N}$, or more generally $x^2 \equiv a \pmod{p^N}$ for a prime number $p$.
 A: Can you solve it in $p$-adics?
For any $2$-adic integer $n$ we have the familiar Taylor/binomial-power series for $\sqrt{1+8n}$:
$\sqrt{1+8n}=1+(1/2)(8n)-...$
All omitted terms vanish $\bmod 512$ when $n$ is a multiple of $8$, thus with $n=8$
$\sqrt{65}\equiv 1+(1/2)(8×8) \bmod 512$
$\sqrt{65}\equiv 1+32\equiv 33 \bmod 512$
So one square root (the root closer to $1$ in $2$-adics) is $33 \bmod 512$.
For reference here are more terms in the binomial expansion used above:
$\sqrt{1+8n}=1+(2^2)n-(2^3)n^2+(2^5)n^3-(5×2^5)n^4+(7×2^7)n^5-(21×2^8)n^6+(33×2^{10})n^7-(429×2^9)n^8+...$
This is sufficiently many terms to obtain at least eleven bits in the square root ($\bmod 2048$) for any integer $n$ (not just multiples of $8$ as above); the exponent on $2$ is at least $11$ in all subsequent terms.  The exponent on $2$ is always at least the number of terms rendered (e.g. $2^9$ in the ninth term), thus guaranteeing convergence in $2$-adics.  Note also that for a conventional integer square the obtained square root is one greater than a multiple of $4$, which in conventional algebra could be the negative root (e.g. $\sqrt{9}=-3$).
A: A basic refinement would be to use the fact that if $x^2\equiv y \pmod{p^N}$, then for every $N'<N$ we also have $x^2\equiv y\pmod{p^{N'}}$. You might as well assume $y$ is coprime to $p$, since otherwise you can easily pull out its factors of $p$.
You can then proceed as follows: Find every solution $x$ to this mod $p$. Then, each of these solutions could lift to $p$ things mod $p^2$ - so you can find the solutions mod $p^2$ checking those. Then you can lift again to $p^3$ and so on - and, since there will be at most four* square roots mod $p^N$ (and only two if $p$ is odd). Thus, you'll only have to check about $4pN$ values for this method to work - a lot better than checking $p^N$ values - but wait, we can do even better!

Note that the actual lifting step can be further simplified: If we know that $x^2\equiv y\pmod{p^N}$, then we can try looking for some $z=x+p^Nc$ such that $z^2\equiv y\pmod{p^{N+1}}$. Working this out gives $x^2+2p^Nc\equiv y\pmod{p^{N+1}}$, which is easy to solve for $c$ (or, if $p=2$, we find that $c$ doesn't even matter - half of our previous square roots lift to two new square roots, and the other half don't lift at all!). Doing this, once you figure out the square root mod $p$, it's basically clean sailing from there.

As for actually solving $x^2\equiv y\pmod{p}$, there are some ways - for instance, if $p=4k+3$, you can find one using Fermat's little theorem, but it seems you're interesting in small primes (e.g. $2$) where it's not hard to check every possible $x$.

(*I also should explicitly point out that, for $k\geq 3$, there are four solutions to $x^2=1\pmod{2^k}$ - and consequently, any number with a square root actually has four square roots - the ones you missed for $65$ mod $512$ were $\pm 223$, which can be found by adding $256$ to your existing roots, since adding $256$ doesn't change the square mod $512$.)
A: This isn't a general method, but here's a neat trick if $2^k|a-1$ for some large enough $k$.
Claim. If $x^2\equiv a\bmod 2^n$ and $a\equiv 1\bmod 2^k$, with $n\leq 2k$, then
$$x\equiv \pm \frac{a+1}{2}\bmod 2^{n-1}.$$
Proof. As there are only two residue classes $\bmod 2^{n-1}$ that will satisfy this, it suffices to show that these both do. Setting $b=(a+1)/2$ we need to show that
$$b^2\equiv a=2b-1\bmod 2^n.$$
Indeed, this is equivalent to
$$2^n|(b-1)^2 \Leftrightarrow b\equiv 1\bmod 2^{\left\lceil\frac{n}{2}\right\rceil},$$
which is itself equivalent to
$$a\equiv 1\bmod 2^{\left\lceil\frac{n}{2}\right\rceil},$$
which is true as $\left\lceil\frac{n}{2}\right\rceil \leq k$.
So, since we have $x^2\equiv a\bmod 2^9$ and $a\equiv 1\bmod 2^6$, we have that $x\equiv \pm \frac{a+1}{2} = \pm 33\bmod 2^8.$
A similar result holds in general:

If $p$ is an odd prime with $x^2\equiv a\bmod p^n$ and $a\equiv 1\bmod 2^k$ with $n\leq 2k$, then
$$x\equiv \pm \frac{a+1}{2}\bmod p^n.$$

This can be proven in nearly identical fashion.
