The number of solutions of $ax^2+by^2\equiv 1\pmod{p}$ is $ p-\left(\frac{-ab}{p}\right)$

Using Legendre symbol, the number of solutions can be written as $$ \sum_{y=0}^{p-1} \left(1+\left(\frac{a-aby^2}{p}\right)\right),$$ since if $a-aby^2$ is nonzero square, you have to count two solutions in $x$, and if $a-aby^2$ is zero, then you have to count one solution in $x$ (namely 0), and if $a-aby^2$ is non-square, then no solutions.

Then rewrite the summation of second term as $$ \sum_{y=0}^{p-1} \left(\frac{a-aby^2}{p}\right)=\left(\frac{-ab}{p}\right)\sum_{y=0}^{p-1}\left(\frac{y^2+d}{p}\right),$$ where $d=-b^*$ ($b^*$ is the inverse of $b$ modulo $p$)

The summation on the right, can be rewritten as $$\sum_{y=0}^{p-1} \left(1+\left(\frac{y}{p}\right)\right)\left(\frac{y+d}{p}\right)$$

Then the sum over single Legendre symbol is zero, but it remains to compute $$\sum_{y=0}^{p-1}\left(\frac{y}{p}\right)\left(\frac{y+d}{p}\right)$$

This one involves a trick: using $\left(\frac{y^*}{p}\right)$ instead of $\left(\frac{y}{p}\right)$

Then you will find that $$\sum_{y=0}^{p-1}\left(\frac{y}{p}\right)\left(\frac{y+d}{p}\right)=-1$$

Now, putting all together, you get the desired answer.


Let $k = -b/a$ and $c = 1/a$ (both are nonzero), so that you can rewrite the equation as $x^2-ky^2 = c$.

Call $N(k,c)$ the number of solutions $(x,y)$ to $x^2-ky^2 = c$. It is easy to check that this number only depends on $c$ and wether $k$ is a square, a nonsquare, or $0$, and that $N(0,c) = p(1 + \binom c p)$

Call $M(y,c)$ the number of solutions $(x,k)$ to $x^2-ky^2 = c$. If $y=0$ then again, $M(0,c) = N(0,c)$. And if $y \neq 0$, we clearly have $k = (x^2-c)/y^2$, hence $M(y,c) = p$ (one value of $k$ for each value of $x$)

The total number $T(c)$ of triplets $(x,y,k)$ satisfying the equation is $T(c) = \sum N(k,c) = \sum M(y,c) = M(0,c) + (p-1)p$, hence $\sum_{k \neq 0} N(k,c) = p^2-p$.


Now, if $k$ is a square $u^2$ you can factor the left-hand side to get $(x-uy)(x+uy) = c$. Since $p>2$, the change of variable $(x,y) \to (x-uy,x+uy) = (s,t)$ is invertible, and you get $st = c$, of which there are $p-1$ solutions for $c \neq 0$, (and $2p-1$ for $c=0$).

This leaves you, when $c \neq 0$, with $(p^2-p)-(p-1)^2/2 = (p-1)(p+1)/2$ triplets with $k$ nonsquare, hence $p+1$ solutions for each nonsquare $k$.
And when $c = 0$, with $(p^2-p)-(p-1)(2p-1)/2 = (p-1)/2$ triplets, hence $1$ solution (the trivial $x=y=0$ one) for each nonsquare $k$.