Efficient way to find squares mod a prime power?

Assume we are given the problem of say finding all squares modulo $3^4$. Is there any efficient way to compute this without having to check a ton of cases? For just a prime we can use quadratic reciprocity, but that doesn't do any good here.

A problem like this was asked on an oral exam in my program (the number was $4000$, so one had to apply the Chinese remainder theorem first and then solve two questions like this), so I guess the professor thought that this is something that one should be able to work out on the black board pretty quickly.

The only approach I know is to first find all squares modulo say $9$, which would give $0,1,4,7$, then look at the numbers:

$0,0+9,0+18,1,1+9,1+18,4,4+9,4+18,...$

and figure out which of these liftings are squares mod $27$. This quickly gets extremely annoying. Especially when having to lift these squares to $3^4$.

This solution depends on the structure of the group $U_{p^\ell}$ of units in the ring $\mathbf{Z}_{p^\ell}$. If you have not covered that, then the question was a bit mean IMHO.

If $p>2$ is a prime, then you can apply the following procedure to find all the squares mod $p^\ell$. Let us first find those square that are coprime to $p$. The group $U_{p^\ell}$ is cyclic of order $\phi(p^\ell)=(p-1)p^{\ell-1}$. In a cyclic group of even order exactly one half of the elements are squares. But for an integer $m$ to be a quadratic residue modulo $p^\ell$ it is necessary for $m$ to be a quadratic residue modulo $p$. This already prevents one half of the elements of $U_{p^\ell}$ from being squares. Therefore by the above observation this necessary criterion is also sufficient: $m$ is a square modulo $p^\ell$, if and only if $m$ is a square modulo $p$.

If $p\mid m$, then $m\equiv a^2\pmod{p^\ell}$ implies that $p\mid a$. Let's write $a=a'p$. Then $m\equiv a'^2p^2\pmod{p^\ell}$, and we see that we must have $p^2\mid m$. Therefore $m=p^2m'$, and $m'\equiv a'^2\pmod {p^{\ell-2}}$. This means that finding squares divisible by $p$ in the ring $\mathbf{Z}_{p^\ell}$ is equivalent to finding all the squares in $\mathbf{Z}_{p^{\ell-2}}$. Rinse. Repeat.

When $p=2$ a similar general approach works. This time $U_{2^\ell}$ is a direct product of two cyclic groups. One of order $2$ and the other of order $2^{\ell-2}$. Therefore exactly one quarter of elements of $U_{2^\ell}$ are squares, and they consist of the numbers $m\equiv 1\pmod 8$, because we see this by a direct calculation in the case $\ell=3$, and for $\ell >3$ the argument above can be repeated. Finding the even squares is done in the same way as in the case $p>2$.

The dead-simple way to list the squares mod $3^4 = 81$ is this:

$\pmatrix{0 & 0 + 1 = 1 & 1 + 3 = 4 & 4 + 5 = 9 & 9 + 7 = 16 & 16 + 9 = 25 & 25 + 11 = 36 & 36 + 13 = 49 & 49 + 15 = 64 \cr 64 + 17 = 0 & 0 + 19 = 19 & 19 + 21 = 40 & 40 + 23 = 63 & 63 + 25 = 7 & 7 + 27 = 34 & 34 + 29 = 63 & 63 + 31 = 13 & 13 + 33 = 46 \cr 46 + 35 = 0 & 0 + 37 = 37 & 37 + 39 = 76 & 76 + 41 = 36 & 36 + 43 = 79 & 79 + 45 = 43 & 43 + 47 = 9 & 9 + 49 = 58 & 58 + 51 = 28 \cr 28 + 53 = 0 & 0 + 55 = 55 & 55 + 57 = 31 & 31 + 59 = 9 & 9 + 61 = 70 & 70 + 63 = 52 & 52 + 65 = 36 & 36 + 67 = 22 & 22 + 69 = 10\cr 10 + 71 = 0 & 0 + 73 = 73 & 73 + 75 = 67 & 67 + 77 = 63 & 63 + 79 = 61\cr}$

... and we can stop there, because $(81-x)^2 \equiv x^2 \mod 81$.