【Analysis】First, if $p \equiv 1(\bmod 4)$, then by Euler's criterion, we know $\left(\frac{-1}{p}\right)=1$, where $\left(\frac{d}{p}\right)$ is the Legendre symbol, defined as
$$
\left(\frac{d}{p}\right)=\left\{\begin{array}{ll}
1, & d \text { is a quadratic residue modulo } p; \\
-1, & d \text { is a non-quadratic residue modulo } p; \\
0, & p \mid d .
\end{array}\right.
$$

Thus, there exists $u$ such that $u^{2} \equiv-1(\bmod p)$.
Hence, $x^{2}+y^{2} \equiv x^{2}-(u y)^{2}(\bmod p)$
$$
=(x-u y)(x+u y) \text {. }
$$

Also, when $p \nmid a$, for any $i \in\{1,2, \cdots, p-1\}$, there exists a unique $k_{i}$ such that $i k_{i} \equiv a(\bmod p)$, therefore, each pair $(x, y)$ can be written as
$$
\begin{array}{l}
\left\{\begin{array}{l}
x-u y \equiv i(\bmod p), \\
x+u y \equiv k_{i}(\bmod p)
\end{array}\right. \\
\Rightarrow\left\{\begin{array}{l}
x \equiv \frac{i+k_{i}}{2}(\bmod p), \\
y \equiv \frac{k_{i}-i}{2} u^{-1}(\bmod p)
\end{array}\right.
\end{array}
$$

in the form, with a total of $p-1$ pairs.
When $a=0$,
$x^{2}+y^{2}=0(\bmod p)$
$\Rightarrow x^{2} \equiv(u y)^{2}(\bmod p)$
$\Leftrightarrow x \equiv \pm u y(\bmod p)$.
Thus, for each $x(x \neq 0), y$ has 2 solutions,
$x$ has $p-1$ choices; and when $x \equiv 0(\bmod p)$,
$y \equiv 0(\bmod p)$,

Therefore, there are $2 p-1$ pairs in total.
Of course, $p \equiv-1(\bmod 4)$ is a more complicated case, representing $\left(\frac{-1}{p}\right)=-1$, i.e., -1 is a non-quadratic residue.
Thus, when $a \neq 0$,
(1) For $\left(\frac{a}{p}\right)=1$, there exists $t$ such that
$t^{2} \equiv a(\bmod p)$.
Hence, $x^{2}+y^{2} \equiv a(\bmod p)$
$$
\Leftrightarrow\left(\frac{x}{t}\right)^{2}+\left(\frac{y}{t}\right)^{2}=1(\bmod p)
$$
$\Rightarrow x^{2}+y^{2} \equiv a(\bmod p)$ has the same number of solutions as $x^{2}+y^{2} \equiv$ $1(\bmod p)$.
(2) For $\left(\frac{a}{p}\right)=-1$, by
$\left(\frac{-1}{p}\right)=-1(\bmod p)$
$$
\Rightarrow\left(\frac{-a}{p}\right) \equiv\left(\frac{-1}{p}\right)\left(\frac{a}{p}\right) \equiv 1(\bmod p)
$$
$\Leftrightarrow$ there exists $t$ such that $t^{2}=-a(\bmod p)$.
By $x^{2}+y^{2} \equiv a(\bmod p)$
$$
\Rightarrow\left(\frac{x}{t}\right)^{2}+\left(\frac{y}{t}\right)^{2} \equiv-1(\bmod p)
$$
$\Rightarrow x^{2}+y^{2} \equiv a(\bmod p)$ has the same number of solutions as $x^{2}+y^{2} \equiv$ $-1(\bmod p)$.

Next, we find the number of solutions to $x^{2}+y^{2} \equiv a^{2}+b^{2}(\bmod p)$.

For $(x-a)(x+a) \equiv(b-y)(b+y)(\bmod p)$, there are $(p-1)(p-2)+p-1$ pairs $(b, y)\left(b^{2} \neq y^{2}(\bmod p)\right)$, each corresponding to $p-1$ pairs of $(x, a)$; there are $2 p-1$ pairs $(b, y)\left(b^{2} \equiv y^{2}(\bmod p)\right)$, corresponding to $2(p-1)+1=2 p-1$ pairs of $(x, a)$.
Therefore, the total number of solutions is
$$
\begin{array}{l}
(p-1)^{3}+(2 p-1)^{2} \\
=p^{3}-3 p^{2}+3 p-1+4 p^{2}-4 p+1 \\
=p^{3}+p^{2}-p .
\end{array}
$$

Let $x^{2}+y^{2} \equiv a(\bmod p)$ have $m$ pairs of solutions, where,
$\left(\frac{a}{p}\right)=1$;
$x^{2}+y^{2} \equiv a(\bmod p)$ have $n$ pairs of solutions, where,
$$
\left(\frac{a}{p}\right)=-1 \text {. }
$$

If $(x, y)$ is not congruent to $(0,0)$ modulo $p$, then
$x^{2}+y^{2} \equiv 0(\bmod p)$
$\Leftrightarrow x^{2} \equiv-y^{2}(\bmod p)$
$\Leftrightarrow\left(x y^{-1}\right)^{2} \equiv-1(\bmod p)$.
This contradicts $\left(\frac{-1}{p}\right)=-1$.
Thus, the solution to $x^{2}+y^{2} \equiv 0(\bmod p)$ is only one pair, which is $(0,0)$.
Then it is easy to know
$$
\left\{\begin{array}{l}
\frac{p-1}{2} m^{2}+\frac{p-1}{2} n^{2}+1=p^{3}+p^{2}-p, \\
\frac{p-1}{2}(m+n)+1=p^{2} .
\end{array}\right.
$$

Solving gives $m=p+1, n=p+1$.
Therefore, $x^{2}+y^{2} \equiv a(\bmod p)$ has $p+1$ pairs of solutions.
$$
\begin{array}{l}
\quad \mid\{(x, y) \mid x^{2}+y^{2} \equiv a(\bmod p), x, y \in\{0, \\
1, \cdots, p-1\}\} \mid
\end{array}
$$