Count number of roots of polynomial modulo prime power I found this problem in a number theory course, I am assuming (but not sure) it is supposed to be an application of Hensel's lemma.

For every $n \in \mathbb{N}_0$, determine the number of solutions over $\mathbb{Z}/5^n\mathbb{Z}$ of the equation $y^2 = x^3 + x^2 + 5$.

In the base case $n = 1$, there are four solutions $(0, 0), (3, 1), (3, 4), (4, 0)$. The problem with Hensel's lemma in the multivariate case seems to be, if you lift your root to a higher prime power, you can no longer be sure this new root is unique. So instead I figured for every root, I'd fix each variable separately, and consider the univariate problem in the other respective variable. From this I derived the following results:


*

*$(0, 0)$ doesn't have higher power solutions by fixing either variable.

*$(3, 1)$ has a unique higher power solution by fixing either variable.

*$(3, 4)$ has a unique higher power solution by fixing either variable.

*$(4, 0)$ has a unique higher power solution by fixing the second variable.


However I have no clue how to further deduce the amount of solutions apart from these implicit ones. Just to get an idea of where I was headed I computed all solutions in the case $n = 2$. This way I did notice there are patterns in the solutions:


*

*$(0, 0)$ mod $5$, there are no corresponding solutions mod $25$.

*$(3, 1)$ mod $5$, the solutions are $(3 + 5k, 21 - 5k)$ mod $25$ for each $k$.

*$(3, 4)$ mod $5$, the solutions are $(3 + 5k, 4 + 5k)$ mod $25$ for each $k$.

*$(4, 0)$ mod $5$, the solutions are $(19, 5k)$ mod $25$ for each $k$.


I can sense there must be a connection between the results I got in the univariate case, and the results I have to get in the multivariate case. But I can't seem to figure out how to formally argue, based on the results of Hensel's lemma, why those patterns in the solutions of $n = 2$ turned out the way they are. Neither do I have a clue how to generalize my results towards an arbitrary power of $5$.
 A: COMMENT.- In general, every solution $x_0$ of $f (x) = 0 \mod p$, except when $\color{red}{f '(x_0) \ne 0 \mod p}$, gives a solution of $f (x) = 0 \mod p ^ n$solving successively the equation modulo $p^2,p^3,\cdots p^n$.
Because of $f (x_0) = 0\mod p ^ n\Rightarrow f (x_0) = 0 \mod p$  putting $x=x_0+px_1$ we can solve $$f(x_0+px_1)=0\mod p^2$$ getting after simplifications$$f(x_0)+px_1f'(x_0)=0\mod p^2$$ which implies  $$f(x_0)+px_1f'(x_0)=0\mod p$$ so we get $$x_1=\frac{-f(x_0)}{pf'(x_0)}\mod p$$    For example the given solution $(3,1)$ modulo $5$ of $y^2=x^3+x^2+5$ gives the univariable equation $$1=x^3+x^2+5$$ which is the equation $$f(x)=x^3+x^2+4=0$$ and with the exposed methode one get $x_1=\dfrac{-f(3)}{5f'(3)}=\dfrac{-(27+9+4)}{5(3\cdot3^2+2\cdot3)}=-\dfrac 33=4$.
Thus $(3,4)$ is a solution of $y^2=x^3+x^2+5$ modulo $25$.
In the proposed equation we have to go first to $3$ distinct equations because if $a$ is a square modulo $25$ we must have $f(x)=x^3+x^2+5-a=0\mod 5$ and because of $x^3+x^2\in\{1,0,2\}$ one has 
$$\begin{cases}f_1(x)=x^3+x^2+5=0\\f_2(x)=x^3+x^2+4=0\\f_3(x)=x^3+x^2+3=0\end{cases}$$
Some work gives then the following solutions modulo 25:
$(3,\pm4),(8,\pm9),(13,\pm11),(18,\pm6),(19,0),(23,\pm1)$.
A: I think I figured it out, but I also feel like I'm trying to re-prove results that should probably follow straight out of Hensel's lemma. For completeness' sake, I'll post the version of Hensel's lemma I have here.

Let $p$ be prime, $k \in \mathbb{Z}_{>0}$, $f \in \mathbb{Z}_p[x_1, \dots, x_k], a = (a_1, \dots, a_k) \in \mathbb{Z}_p^k$ en $e \in \mathbb{Z}_{\ge 0}$. Assume
\begin{equation}
f(a) \equiv 0 \bmod p^{2e + 1}, \quad \min_{1 \le i \le k} \text{ord}_p\left(\frac{\partial f}{\partial x_i}(a)\right) = e.
\end{equation}
  Then some $b \in \mathbb{Z}_p^k$ exists such that $f(b) = 0$ and $b \equiv a \bmod p^{e+1}$.

And this is the solution I came up with.
First, consider the solutions for $k = 1$. Since $y^2$ is a square, it follows that $y \in \{0, 1, 4\}$. There are exactly four solutions $(x, y) \in (\mathbb{Z}/{5\mathbb{Z}})^2$ which are given by $z_1 = (0, 0), z_2 = (4, 0), z_3 = (3, 1), z_4 = (3, 4)$. Now consider the preconditions of the multivariate lemma of Hensel-Rychlik on the polynomial given by $f(x, y) = x^3 + x^2 + 5 - y^2$.
\begin{array}{c|c|c|c|c|c}
   z     & (x, y) & (0, 0) & (4, 0) & (3, 1) & (3, 4) \\
   \hline
   \frac{\delta f}{\delta x}(z) & x(3x + 2) & 0  & 1  & 3  & 3 \\
   \frac{\delta f}{\delta y}(z) & -2y  & 0  & 0  & 3  & 2 \\
   \hline
   \text{min ord}   & -   & \infty & 0  & 0  & 0
  \end{array}
Every zero except the first satisfies the preconditions, which implies at least one solution exists congruent with every zero, except (not necessarily) for $(0, 0)$. Now assume a zero $(5u, 5v) \in \mathbb{Z}/{5\mathbb{Z}}$ of the polynomial exists. Then it follows that
\begin{align*}
 f(5u, 5v) &\equiv (5u)^3 + (5u)^2 + 5 - (5v)^2 \bmod 25 \\
   &\equiv 5 \bmod 25 \equiv 0 \bmod 25
\end{align*}
This leads to a contradiction, which means a solution congruent to $(0, 0)$ is indeed impossible. Now assume some arbitrary $a, b \in \mathbb{Z}/{5^{k + 1}\mathbb{Z}}$ for which holds that $f(a, b) \equiv 0 \bmod 5^{k + 1}$. In other words, it is congruent to one of the previously found roots. Then for some $u, v \in \mathbb{Z}/{5^{k + 1}\mathbb{Z}}$ it follows that
\begin{align*}
 f(a + 5^ku, b + 5^kv) &\equiv (a + 5^ku)^3 + (a + 5^kv)^2 + 5 - (b + 5^kv)^2 \bmod 5^{n + 1} \\
     &\equiv a^3 + 3 \cdot a^25^ku + a^2 + 2 \cdot a5^ku + 5 - b^2 - 2 \cdot b5^kv \bmod 5^{n + 1} \\
     &\equiv 3 \cdot a^25^ku + 2 \cdot a5^ku - 2 \cdot b5^kv \bmod 5^{n + 1} \\
     &\equiv 5^k(3 \cdot a^2u + 2 \cdot au - 2 \cdot bv) \bmod 5^{n + 1} \\
     &\equiv 0 \bmod 5^{n + 1} \\
     &\,\Updownarrow \\
 (3a^2 + 2a)u  &\equiv 2bv \bmod 5
\end{align*}
In other words, the existence of more solutions in $\mathbb{Z}/{5^{k + 1}\mathbb{Z}}$ completely depends on their behavior in $\mathbb{Z}/{5^k\mathbb{Z}}$, which can be checked using the original three roots.


*

*(a, b) = (4, 0): This implies $u \equiv 0 \bmod 5$, which means 5 solutions $(u, v)$ exist.

*(a, b) = (3, 1): This implies $3u \equiv 2v \bmod 5$, which means 5 solutions $(u, v)$ exist.

*(a, b) = (3, 4): This implies $3u \equiv 3v \bmod 5$, which means 5 solutions $(u, v)$ exist.


This implies, for every solution of the equation in $\mathbb{Z}/{5^k\mathbb{Z}}$, five congruent solutions exist in $\mathbb{Z}/{5^{k + 1}\mathbb{Z}}$. It can be concluded that, for $k > 1$, the amount of solutions to the equation is given by $3 \cdot 5^k$.
