Note: Definitions for certain terms are posted below.
It is quite easy if you use the ring $\mathbb{Z}[\text{i}]$ of Gaussian integers. For $w\in\mathbb{Z}[\text{i}]$, write $\bar{w}$ for the complex conjugate of $w$. Factorize $z\in \mathbb{N}$ as a product of primes $$z=2^mp_1^{r_1}p_2^{r_2}\cdots p_k^{r_k}q_1^{s_1}q_2^{s_2}\cdots q_l^{s_l}$$ over $\mathbb{Z}$, where the $p_i$'s and the $q_j$'s are positive odd primes such that $p_1,p_2,\ldots,p_k\equiv 1\pmod{4}$ and $q_1,q_2,\ldots,q_l\equiv -1\pmod{4}$. It is well known that $z$ is a sum of two square if and only if each of the $s_j$'s is even. From now on, we assume that $s_j$ is even for every $j=1,2,\ldots,l$.
For $i=1,2,\ldots,k$, $p_i$ can be factorized into the product $\pi_i\bar{\pi}_i$, where $\pi_i\in\mathbb{Z}[\text{i}]$ is a Gaussian prime. Also, $2=(1+\text{i})(1-\text{i})$. Hence, you can pick $\mu_i\in\left\{0,1,\ldots,r_i\right\}$ for each $i=1,2,\ldots,k$ and write $$z=\left((1+\text{i})^m\,\prod_{i=1}^k \pi_i^{\mu_i}\bar{\pi}_i^{r_i-\mu_i}\,\prod_{j=1}^l q_j^{s_j/2}\right)\left((1-\text{i})^m\,\prod_{i=1}^k \bar{\pi}_i^{\mu_i}\pi_i^{r_i-\mu_i}\,\prod_{j=1}^l q_j^{s_j/2}\right)\,.$$ Hence, if $x,y\in\mathbb{Z}$ are such that $$x+y\text{i}=u\,(1+\text{i})^m\,\prod_{i=1}^k \pi_i^{\mu_i}\bar{\pi}_i^{r_i-\mu_i}\,\prod_{j=1}^l q_j^{s_j/2}\,,$$ where $u\in\{-1,+1,-\text{i},+\text{i}\}$, then $x-y\text{i}=\bar{u}\,(1-\text{i})^m\,\prod_{i=1}^k \bar{\pi}_i^{\mu_i}{\pi}_i^{r_i-\mu_i}\,\prod_{j=1}^l q_j^{s_j/2}$ and
$$z=(x+y\text{i})(x-y\text{i})=x^2+y^2\,.$$
It can be shown that there are $4\prod_{i=1}^k\left(r_i+1\right)$ ways to find $(x,y)\in\mathbb{Z}^2$ such that $x^2+y^2=z$. Since there are $4$ ways to choose $u$ and, for $i=1,2,\ldots,k$, there are $r_i+1$ ways to choose $\mu_i$, we have found all possible pairs $(x,y)\in\mathbb{Z}^2$ such that $x^2+y^2=z$. If signs and ordering of $x,y$ are ignored, there are exactly $\left\lceil\frac{1}{2}\prod_{i=1}^k\left(r_i+1\right)\right\rceil$ possible pairs. If you want to exclude the case where $x$ or $y$ may be zero, then there are two cases:
(1) All $r_i$'s are even: there are $\frac12\left(\prod_{i=1}^k\left(r_i+1\right)-(-1)^m\right)$ pairs of nonzero $x,y$, up to signs and order;
(2) At least one of the $r_i$'s is odd: there are $\frac{1}{2}\prod_{i=1}^k\left(r_i+1\right)$ pairs of nonzero $x,y$, up to signs and order.
For example, let $z:=90=2\cdot 5\cdot 3^2$, with $2=(1+\text{i})(1-\text{i})$ and $5=(2+\text{i})(2-\text{i})$. Hence, we can take $x,y\in\mathbb{Z}$ such that $x+y\text{i}=(1+\text{i})\cdot (2+\text{i})\cdot 3=3+9\text{i}$, or $(x,y)=(3,9)$. In fact, there are $8$ ways to write $z=90$ as a sum of two squares; i.e., $(\pm3,\pm9)$ and $(\pm9,\pm3)$ are all possible values of $(x,y)$ (only one pair---$(x,y)=(3,9)$---if sings and ordering are ignored). Another example is $z:=50$. There are $12$ ways to write $z$ as a sum of two squares; i.e., $(\pm5,\pm5)$, $(\pm1,\pm7)$, and $(\pm7,\pm1)$ are all possible values of $(x,y)$ (only two pairs---$(x,y)=(5,5)$ and $(x,y)=(1,7)$---if signs and ordering are ignored).
In the case where $z=t^2$ for some $t\in\mathbb{N}$, the situation is not changed at all. You still need to factorize $z$ as above. For example, if $z=13^2\cdot 17^2\cdot 19^2$, you notice that $13=(3+2\text{i})(3-2\text{i})$ and $17=(4+\text{i})(4-\text{i})$. So, you can take $x,y\in\mathbb{Z}$ such that $x+y\text{i}=(3+2\text{i})\cdot(4+\text{i})\cdot 19=19(10+11\text{i})$, or $(x,y)=(10\cdot 19,11\cdot 19)$. Indeed, there are $36$ possible solutions $(x,y)\in\mathbb{Z}^2$ to $x^2+y^2=z$, and, up to signs and order, there are $5$ of them: $(10\cdot 19,11\cdot 19)$, $(5\cdot 19,14\cdot 19)$, $(1\cdot 13\cdot 19,4\cdot 13\cdot 19)$, $(2\cdot 17\cdot 19,3\cdot 17\cdot 19)$, and $(0,13\cdot 17\cdot 19)$.
P.S. If you are asking how to factorize an integer prime $p\in\mathbb{N}$ such that $p\equiv 1\pmod{4}$ into a product $\pi\bar{\pi}$ of conjugate Gaussian primes $\pi$ and $\bar{\pi}$, then you first need to solve for $t^2\equiv -1\pmod{p}$ (you can, in fact, take $t:=\left(\frac{p-1}{2}\right)!$, but it is probably easier to search for $t\in\left\{1,2,\ldots,\frac{p-1}{2}\right\}$ directly, when $p$ is large). The ring $\mathbb{Z}[\text{i}]$ has a division algorithm. Therefore, you can define a Euclidean algorithm in $\mathbb{Z}[\text{i}]$. Then, you take $\pi$ to be a greatest common divisor of $p$ and $t+\text{i}$.
For example, if $p=29$, we take $t:=12 \equiv \left(\frac{29-1}{2}\right)!\pmod{29}$. Perform the Euclidean algorithm to get $p=(t+\text{i})\cdot 2+(5-2\text{i})$ and $t+\text{i}=(5-2\text{i})\cdot(2+\text{i})+0$, from which we conclude that $5-2\text{i}$ is a greatest common divisor of $p$ and $t+\text{i}$. Hence, $\pi$ can be taken to be $5-2\text{i}$. Note that $p=\pi\bar{\pi}$, as desired. For more information, see the notes by Keith Conrad below.
Definitions
Factorize: "resolve or be resolvable into factors." Google
Gaussian integer: "a complex number whose real and imaginary parts are both integers" Wikipedia
Reference
http://www.math.uconn.edu/~kconrad/blurbs/ugradnumthy/Zinotes.pdf
$18^2 + 24^2 = 324 + 576 = 900 = 30^2$
– bthmas Jun 28 '15 at 06:04