How can I count the number of solutions to $x^2 + y^2 \equiv 1 \pmod{p}$? I know that the result is that there are $p+1$ solutions if $p \equiv 3 \pmod{4}$, and $p-1$ solutions if $p \equiv 1 \pmod{4}$, but I have no idea how to prove this. I've written a program to verify that this is the case, but it hasn't helped me prove it. I know there are always the following 4 solutions, $(0,1)$, $(1,0)$, $(p-1,0)$, $(0,p-1)$.
 A: The number of solutions to $x^2\equiv a\pmod p$ is $1+\left(\frac ap\right)$ (Legendre symbol). So the number of solutions of $x^2+y^2\equiv1\pmod{p}$ is
$$S=p+\sum_{x=0}^{p-1}\left(\frac{1-x^2}{p}\right)
=p+\sum_{z=1}^{p-1}\left(\frac{-2z-z^2}{p}\right)
=p+\sum_{z=1}^{p-1}\left(\frac{-2\overline z-1}{p}\right)$$
where I have put $x=1+z$ and $\overline z$ is the inverse of $z$ modulo $p$. As $z$ varies over $\{1,\ldots,p-1\}$ so does $\overline z$.
Therefore $-2\overline z-1$ ranges over all resides modulo $p$
except for $p-1$.
Therefore
$$S=p+\sum_{t=0}^{p-1}\left(\frac tp\right)-\left(\frac{-1}p\right)=p
-\left(\frac{-1}p\right).$$
An alternative argument in to use the fact that the number of points
on a nonsingular conic in the projective plane over $\Bbb F_p$
is $p+1$. The conic $X^2+Y^2=Z^2$ has two points "at infinity"
when $\left(\frac{-1}p\right)=1$ but none when $\left(\frac{-1}p\right)=-1$.
A: The solutions are precisely the points on the circle that have coordinates in $\Bbb{F}_p$. These can be parametrised just like the rational points on the circle; you know that $(0,1)$ is a solution, so the points on the circle are parametrized by the lines through $(0,1)$ with slope $\lambda\in\Bbb{F}_p$. Explicitly the line through $(0,1)$ with slope $\lambda$ is given by
$$L_{\lambda}:=\{(x,1+\lambda x):\ x\in\Bbb{F}_p\},$$
and the point $(x,1+\lambda x)$ is on the circle if and only if
$$x^2+(1+\lambda x)^2\equiv1\pmod{p}.$$
The trivial solution $x=0$ yields the initial point $(0,1)$, and for $x\neq0$ we can reduce this to
$$(1+\lambda^2)x\equiv-2\lambda\pmod{p}.$$
If $p\equiv3\pmod{4}$ then $1+\lambda^2\neq0$ so we get precisely $p$ more solutions of the form $x=-\frac{2\lambda}{1+\lambda^2}$, yielding a total of $p+1$ solutions.
If $p\equiv1\pmod{4}$ then there are two values of $\lambda$ such that $1+\lambda^2=0$, so we get only $p-2$ more solutions of the form $x=-\frac{2\lambda}{1+\lambda^2}$, yielding a total of $p-1$ solutions.
