Since this old question has resurfaced, let me sketch two ways to prove the stated formulæ using algebraic geometry:

The first way is fairly elementary. Let us stick for definiteness with the number $f(n) := \#\{(x,y,z) \in (\mathbb{Z}/n\mathbb{Z})^3 : x^2+y^2+z^2=1\}$ of points on the ($2$-)sphere in $3$-dimensional affine space $\mathbb{A}^3_{\mathbb{Z}/n\mathbb{Z}}$ over $\mathbb{Z}/n\mathbb{Z}$. As others have pointed out, by the Chinese Remainder Theorem, $f$ is multiplicative. Furthermore, as $\{(x,y,z)\in\mathbb{A}^3_{\mathbb{Z}} : x^2+y^2+z^2=1\}$ is smooth outside of the prime $2$ (because the equation and its partial derivatives generate the unit ideal in $\mathbb{Z}[\frac{1}{2}][x,y,z]$), an appropriate form of Hensel's lemma or the definition of (formal) smoothness ensures that $f(p^{k+1}) = p^2\cdot f(p^k)$ whenever $p$ is an odd prime (the exponent $2$ in $p^2$ is the dimension of the tangent space).

So, leaving aside the case of the prime $2$ (I'm merely sketching the argument here), we are left with proving that the number $f(p)$ of points over $\mathbb{F}_p$ is $p(p+1)$ or $p(p-1)$ according as $p$ is congruent to $1$ or $3$ mod $4$. Now stereographic projection $\varphi\colon (x,y,z) \mapsto (\frac{x}{1-z}, \frac{y}{1-z})$ defines a birational map between the sphere and the plane, with inverse $(u,v) \mapsto (\frac{2u}{1+u^2+v^2}, \frac{2v}{1+u^2+v^2}, \frac{-1+u^2+v^2}{1+u^2+v^2})$, so it is "almost" a bijection, and we just have to analyse the points where it fails: $\varphi$ is defined everywhere except when $z=1$ which, when $p\equiv 3\pmod{4}$, is just one point $(0,0,1)$, whereas when $p\equiv 1\pmod{4}$, consists of the $(x,\pm\sqrt{-1}\,x,1)$, so, $2p-1$ of them; as for the image of $\varphi$, it consists of those $(u,v)$ such that $1+u^2+v^2 \neq 0$, which is the complement of a conic having $p-1$ or $p+1$ points according as $p$ is congruent to $1$ or $3$ mod $4$, so the image of $\varphi$ has $p^2-p+1$ or $p^2-p-1$ points, and finally the sphere has $p^2+p$ or $p^2-p$ in each case.

The second, more sophisticated way, works in any dimension. Once again leaving aside the troublesome prime $2$, the idea is to write the ($(n-1)$-)sphere in dimension $n$ as $\mathit{SO}(n)/\mathit{SO}(n-1)$ where $\mathit{SO}(n)$ refers to the orthogonal group of the quadratic form $x_1^2 + \cdots + x_n^2$. Now the latter is a "standard" quadratic form over the reals, but over $\mathbb{F}_p$, when $n=2m$ is even, this form is said to be of "plus type" (= Witt index $m$) for $p\equiv 1\pmod{4}$ or $m$ even, and of "minus type" (= Witt index $m-1$) for $p\equiv 3\pmod{4}$ and $m$ odd; for $n$ odd, there is only one nondegenerate quadratic form over $\mathbb{F}_q$. Now the order of the groups $\mathit{SO}(2m,{+})$ (for a quadratic form of "plus" type), $\mathit{SO}(2m,{-})$ ("minus" type) and $\mathit{SO}(2m+1)$ over $\mathbb{F}_q$ are known (up to trivial factors, they are what group theorists write $D_m$, $^2D_m$ and $B_m$), and are one half of:

$\#\mathit{Spin}(2m,{+})(\mathbb{F}_q) = q^{m(m-1)} (q^m-1) \prod_{i=1}^{m-1}(q^{2i}-1)$

$\#\mathit{Spin}(2m,{-})(\mathbb{F}_q) = q^{m(m-1)} (q^m+1) \prod_{i=1}^{m-1}(q^{2i}-1)$

$\#\mathit{Spin}(2m+1)(\mathbb{F}_q) = q^{m^2} \prod_{i=1}^{m}(q^{2i}-1)$

(See any number of books on finite simple groups or algebraic groups for a proof and explanation of these formulæ.) So by taking the ratio of these quantities we get the number of points on the sphere: for example, for the ($2$-)sphere in dimension $3$, we have $\#\mathit{Spin}(3)(\mathbb{F}_p) = p^3 - p$ and $\#\mathit{Spin}(2,{+})(\mathbb{F}_p) = p - 1$ and $\#\mathit{Spin}(2,{-})(\mathbb{F}_p) = p + 1$, giving $f(p) = p(p+1)$ for $p\equiv 1\pmod{4}$ and $f(p) = p(p-1)$ for $p\equiv 3\pmod{4}$. But of course, the algebraic group method makes it possible to compute many more things.

3more comments