$p \equiv 1 \pmod{4}$. Show that there are exactly $4(k+1)$ integer pairs $(a,b)$ s.t. $a^2+b^2 =p^k$ 
$p \equiv 1 \pmod 4$. Show that there are exactly $4(k+1)$ integer pairs $(a,b)$ s.t. $a^2+b^2 =p^k$

I have proved the less general statement for $k=1$ using the fact that the Gaussian integers are a integral domain and factorising it into irreducible then changing this factorisation multiplying but units. 
This doesn’t help in this case because $p^k$ doesn’t factorise into irreducible so I’m not sure what to do. 
 A: Since $p=\pi\bar{\pi}$ for some Gaussian prime $\pi\in R:=\mathbb{Z}[\text{i}]$ and $$(a+b\text{i})(a-b\text{i})=p^k=\pi^k\bar{\pi}^k\,.\tag{*}$$  By the unique factorization property of $R$, $$a+b\text{i}=u\pi^s\bar{\pi}^{t}$$ for some unit $u$ of $R$ and for some integers $s$ and $t$ such that $0\leq s,t\leq k$.
Now,
$$a-b\text{i}=\overline{a+b\text{i}}=\overline{u\pi^s\bar{\pi}^t}=\bar{u}\bar{\pi}^s\pi^t\,.$$
By (*), we obtain
$$p^{s+t}=1\cdot p^{s+t}=(u\bar{u})\cdot(\pi\bar{\pi})^{s+t}=(a+b\text{i})(a-b\text{i})=p^k\,.$$
That is, $s+t=k$, or $t=k-s$.  Hence,
$$a+b\text{i}=u\pi^s\bar{\pi}^{k-t}$$
with $s\in\{0,1,2,\ldots,k\}$, and $u\in\{+1,-1,+\text{i},-\text{i}\}$.  There are $4$ choices for $u$ and $k+1$ choices for $s$.
In general, suppose that
$$N=2^\alpha\,\left(\prod_{i=1}^m\,p_i^{\beta_i}\right)\,\left(\prod_{j=1}^n\,q_j^{\gamma_j}\right)\,,$$
where $\alpha$, $\beta_1,\beta_2,\ldots,\beta_m$, and $\gamma_1,\gamma_2,\ldots,\gamma_n$ are nonnegative integers, and $p_1,p_2,\ldots,p_m$ are pairwise distinct prime natural numbers congruent to $1$ modulo $4$, and $q_1,q_2,\ldots,q_n$ are pairwise distinct prime natural numbers congruent to $3$ modulo $4$.  If $f(N)$ denotes the number of solutions $(a,b)\in\mathbb{Z}\times\mathbb{Z}$ to the equation $$a^2+b^2=N\,,\tag{$\star$}$$
then
$$f(N)=\left\{
\begin{array}{ll}
0&\text{if }\gamma_j\text{ is odd for some }j\in\{1,2,\ldots,n\}\,,
\\4\,\prod\limits_{i=1}^m\,(\beta_i+1)&\text{if }\gamma_j\text{ is even for every }j\in\{1,2,\ldots,n\}\,.
\end{array}
\right.$$
If $g(N)$ is the number of solutions $(a,b)\in\mathbb{Z}_{\geq 0}\times\mathbb{Z}_{\geq 0}$ to ($\star$), then
$$g(N)=\left\{
\begin{array}{ll}
0&\text{if }\gamma_j\text{ is odd for some }j\in\{1,2,\ldots,n\}\,,
\\\prod\limits_{i=1}^m\,(\beta_i+1)&\text{if }\gamma_j\text{ is even for every }j\text{ and }N\text{ is not a square}\,,\\
\prod\limits_{i=1}^m\,(\beta_i+1)+1&\text{if }\gamma_j\text{ is even for every }j\text{ and }N\text{ is a square}\,.
\end{array}\right.$$
If $g'(N)$ is the number of solutions $(a,b)\in\mathbb{Z}_{> 0}\times\mathbb{Z}_{> 0}$ to ($\star$), then
$$g'(N)=\left\{
\begin{array}{ll}
0&\text{if }\gamma_j\text{ is odd for some }j\in\{1,2,\ldots,n\}\,,
\\\prod\limits_{i=1}^m\,(\beta_i+1)&\text{if }\gamma_j\text{ is even for every }j\text{ and }N\text{ is not a square}\,,\\
\prod\limits_{i=1}^m\,(\beta_i+1)-1&\text{if }\gamma_j\text{ is even for every }j\text{ and }N\text{ is a square}\,.
\end{array}\right.$$
If $h(N)$ is the number of solutions $(a,b)\in\mathbb{Z}_{\geq 0}\times\mathbb{Z}_{\geq 0}$ to ($\star$) such that $a\leq b$, then
$$h(N)=\left\{
\begin{array}{ll}
0&\text{if }\gamma_j\text{ is odd for some }j\in\{1,2,\ldots,n\}\,,
\\
\frac{\prod\limits_{i=1}^m\,(\beta_i+1)}{2}&\text{if }\gamma_j\text{ is even for every }j\text{ and }\beta_i\text{ is odd for some }i\,,
\\
\frac{\prod\limits_{i=1}^m\,(\beta_i+1)+1}{2}&\text{if }\gamma_j\text{ is even for every }j\text{ and }\beta_i\text{ is even for every }i\,.
\end{array}
\right.$$
If $h'(N)$ is the number of solutions $(a,b)\in\mathbb{Z}_{> 0}\times\mathbb{Z}_{> 0}$ to ($\star$) such that $a\leq b$, then
$$h'(N)=\left\{
\begin{array}{ll}
0&\text{if }\gamma_j\text{ is odd for some }j\in\{1,2,\ldots,n\}\,,
\\
\frac{\prod\limits_{i=1}^m\,(\beta_i+1)}{2}&\text{if }\gamma_j\text{ is even for every }j\text{ and }\beta_i\text{ is odd for some }i\,,
\\
\frac{\prod\limits_{i=1}^m\,(\beta_i+1)-(-1)^\alpha}{2}&\text{if }\gamma_j\text{ is even for every }j\text{ and }\beta_i\text{ is even for every }i\,.
\end{array}
\right.$$
A: Let:


*

*$P$ be the set of integer pairs $(a, b)$ with $a^2+b^2 = p^k$;

*$E$ the set of elements in $\mathbb Z[i]$ with norm $p^k$;

*$I$ the set of ideals in $\mathbb Z[i]$ with norm $p^k$.


Because $\mathbb Z[i]$ has unique factorization of ideals, $|I| = k+1$.
We have maps $f : P \to E : (a, b) \mapsto a+bi$ and $g : E \to I : \alpha \mapsto (\alpha)$.
It is clear that $f$ is a bijection. Because $\mathbb Z[i]$ is a PID, $g$ is surjective. The fibers of $g$ are exactly the orbits under multiplication by the unit group, that is, $\{1, i, -1, -i\}$. Thus $|P| = |E| = 4 |I|$.
