counting points on unit sphere mod p Let $f(n)$ be the number of points on the unit sphere $x^2 + y^2 + z^2 = 1\; \pmod n$ with $x,y,z \in \mathbb{Z}/n\mathbb{Z}$  
This is sequence A087784 in the Online Encyclopedia of Integer sequences:

1, 4, 6, 24, 30, 24, 42, 96, 54, 120...

There is a (due to Bjorn Poonen) indicating some regularity to the solutions to this congurence
$$f(n) = n^2* \left\{\begin{array}{cl}3/2&\text{if}\quad\quad 4|n \\
1 &\text{otherwise}
\end{array}\right\}*\prod_{\substack{p|n \\ 1 \mod 4}} \left( 1 + \frac{1}{p}\right)* \prod_{\substack{p|n \\ 3 \mod 4}} \left( 1 - \frac{1}{p}\right)$$
What are some proofs to this identity ?  

Sequence A060968 is the number of points on the unit circle $x^2 + y^2  \equiv 1\; \pmod n$ 

1, 2, 4, 8, 4, 8, 8, 16, 12, 8, 12,...

with a similar multiplicative formula:
 $$g(n) = n* \left\{\begin{array}{cl}2&\text{if}\quad\quad 4|n \\
1 &\text{otherwise}
\end{array}\right\}*\prod_{\substack{p|n \\ 1 \mod 4}} \left( 1 + \frac{1}{p}\right)* \prod_{\substack{p|n \\ 3 \mod 4}} \left( 1 - \frac{1}{p}\right)$$
Perhaps there is a tower of such identities.

The multiplicative structure of these formulas could have an algorithmic interpretation.  The formula for the Euler phi function
\[ \phi(n) =  n \prod_{p|n} \left( 1 - \frac{1}{p} \right) \]
This suggests a sieving algorithm to generate the list of numbers relatively prime to n


*

*write down the numbers $\{ 1, 2, \dots, n \}$

*for reach prime $p|n$ cross out multiples


I'd be especially interested if this type of algorithm existed for $f(n), g(n)$.
 A: 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.
A: It was mentioned by Abhinav Kumar that $f(n)$ is a multiplicative function. So it is sufficient to prove desired formula for $n=p^\alpha$ where $p$ is a prime. Let $p\ge 3$. Then
\begin{gather*}
f(n)=\frac{1}{ p^\alpha}\sum_{t=1}^{p^\alpha}\sum_{x,y,z=1}^{p^\alpha}e_{p^\alpha}((x^2+y^2+z^2-1)t),
\end{gather*}
where $e_n(x)=e^{2\pi i x/n}$.
Suppose that $t=p^\beta\tau$ with $(\tau,p)=1$ and $0\le \beta\le \alpha$. Then
\begin{gather*}
f(n)=\frac{1}{ p^\alpha}\sum_{\beta=0}^{\alpha}\sum_{\tau=1,(\tau,p)=1}^{p^{\alpha-\beta}}\sum_{x,y,z=1}^{p^\alpha}
e_{p^{\alpha-\beta}}((x^2+y^2+z^2-1)\tau)=\\=
\frac{1}{ p^\alpha}\sum_{\beta=0}^{\alpha}p^{3\beta}\sum_{\tau=1,(\tau,p)=1}^{p^{\alpha-\beta}}\sum_{x,y,z=1}^{p^{\alpha-\beta}}
e_{p^{\alpha-\beta}}((x^2+y^2+z^2-1)\tau)=\\=\frac{1}{ p^\alpha}\sum_{\beta=0}^{\alpha}p^{3\beta}\sum_{\tau=1,(\tau,p)=1}^{p^{\alpha-\beta}}
e_{p^{\alpha-\beta}}(-\tau)S^3(\tau,p^{\alpha-\beta}),
\end{gather*}
where $S(a,n)$ is a Gauss sum
$$S(a,n)=\sum_{x=1}^{n}e_n(ax^2).$$ It is known that for odd $n$
$$S(a,n)=\varepsilon_n\left(\frac{a}{ n}\right)\sqrt{n},$$
where $\varepsilon_n=1$ for $n\equiv 1\pmod{ 4}$ and $\varepsilon_n=i$ for $n\equiv 3\pmod{ 4}$
(see Korobov N. M. Exponential sums and their applications, 1992). It means that
\begin{gather*}
f(n)=p^{\frac{\alpha}{2 }}\sum_{\beta=0}^{\alpha}p^{\frac{3\beta}{2 }}\varepsilon^3_{p^{\alpha-\beta}}
\sum_{\tau=1,(\tau,p)=1}^{p^{\alpha-\beta}}
e_{p^{\alpha-\beta}}(-\tau)\left(\frac{\tau}{ p}\right)^{\alpha-\beta}.
\end{gather*}
In the case $\beta=\alpha$ we'll get $p^{2\alpha}$. In the case $\beta=\alpha-1$ we'll have the sum
\begin{gather*}
p^{\frac{\alpha}{2 }+\frac{3(\alpha-1)}{2 }}\varepsilon_p^3\sum_{\tau=1}^{p-1}e_p(-\tau)\left(\frac{\tau}{ p}\right)=
\\=p^{2\alpha-\frac{3}{2 }}\varepsilon_p^3S(-1,p)=p^{2\alpha-1}\left(\frac{-1}{ p}\right).
\end{gather*}
If $\beta<\alpha-1$ then the inner sum vanishes. Let $\tau=\tau_1+p\tau_2$, $1\le \tau_1<p$, $1\le \tau_2\le p^{\alpha-\beta-1}$. We have
\begin{gather*}
\sum_{\tau=1,(\tau,p)=1}^{p^{\alpha-\beta}}
e_{p^{\alpha-\beta}}(-\tau)\left(\frac{\tau}{ p}\right)^{\alpha-\beta}=\sum_{\tau_1=1}^{p-1}\left(\frac{\tau_1}{ p}\right)^{\alpha-\beta}\sum_{\tau_2=1}^{p^{\alpha-\beta-1}}
e_{p^{\alpha-\beta}}(-(\tau_1+p\tau_2))=
\sum_{\tau_1=1}^{p-1}\left(\frac{\tau_1}{ p}\right)^{\alpha-\beta}
e_{p^{\alpha-\beta}}(-\tau_1)\sum_{\tau_2=1}^{p^{\alpha-\beta-1}}
p^{\alpha-\beta-1}\delta_{p^{\alpha-\beta-1}}(1)=0,
\end{gather*}
where $\delta_n(x)$ is the characteristic function of divisibility by $n$:  $\delta_n(x)=1$ for $x\equiv 0\pmod{n }$ and $\delta_n(x)=0$ for $x\not\equiv 0\pmod{n }$.
For $p=2$ one should apply similar formula for Gauss su: $S(a,n)=0$ if $n\equiv 2\pmod{ 4}$ and
$$S(a,n)=(1+i)\varepsilon_n^{-1}\left(\frac{n}{a}\right)\sqrt{n},$$ see Generalized quadratic Gauss sums.
