distribution of $\sqrt{-1} \mod p$ While reading up on quadratic reciprocity, I learned that if $p = 4k+1$ then $-1$ has a square root in $\mathbb{Z} / p \mathbb{Z}$.
Let $r_p$ be an integer with $0\leq r_p < p$ and $r_p^2 \equiv -1 \mod p$. How then is $\frac{r_p}{p} \in \mathbb{Q}$ distributed in $[0,1]$?  Naively I would guess this is uniform distribution.  How can we prove that?

Edit I noticed in the comments, it might be simpler to ask about the equidistribution of $$\{ \tfrac{1}{\sqrt{p}}(a,b): a^2 + b^2 = p\} \subset S^1$$
still in the case $p = 4k+1$.
 A: The equidistribution of the roots of quadratic congruences $\pmod p$ (such as $x^2+1$ in the question) was established in a famous paper of Duke, Friedlander and Iwaniec.  The proof uses sieve ideas as well as ideas from the theory of modular forms. 
A: There is a related result of Friedlander-Iwaniec embedded in the proof of lemma 3.2, which is not equidistribution, but gives an extremely simple heuristic reason why these fractions are well-distributed. Namely, if we fix $X$, restrict to an interval of integers $8X/9 < d \le X$, and consider the set
$S=\{\nu/d : \nu ^2 +1 \equiv 0 \pmod d, d\in(8X/9,X]\}$,
then these points repel each other $\pmod 1$: while we would expect their differences to be on the order of $1/X^2$, it turns out that every pair is at least $1/4X$ apart. 
The idea is that if $\nu^2 +1 \equiv 0\pmod d$ then $\nu=rs^{-1}\pmod d$ where $r,s$ satisfy $r^2+s^2 = d$ exactly. But then $\nu/d\pmod{1}$ is very close to $-\bar{r}/s$, where $\bar{r}$ is the inverse $\pmod s$ of $r$. Now these fractions $-\bar{r}/s$ all have much smaller denominators $<\sqrt{X}$, and so must be on the order of $1/X$ apart.
