How do solutions to this quadratic congruence distribute as the number of factors grows? This is a question that arose during a conversation with a colleague regarding Landau's fourth problem, which asks whether there are infinitely many primes of the form $n^2+1$. The Conjectured asymptotic for
$$\sum_{n\leq X}\Lambda(n^2+1)$$
is 
$$-\sum_{n\leq x^2+1}\frac{\mu(n)\rho(n)\log n}{n}\sim \frac{CX}{\log X}$$
where $\rho(n)$ counts the number of solutions to the congruence $a^2+1\equiv 0 \pmod n$ and $C$ is specific arithmetic constant. Although Landau's fourth problem is considered to be out of reach, the question I want to ask seems to be a natural prerequisite. I have not proved this implication yet, but it seems to be a nice problem anyway. 
Question
 Let $n\in\mathbb{N}$ be squarefree and such that $p|n\Rightarrow p=4m+1$ for some $m\in\mathbb{N}$ so $-1$ is a square $\pmod n$ with $2^{\omega(n)}$ square roots $\pmod n$. Do these roots equidistribute $\pmod n$ as $\omega(n)\rightarrow\infty$?
Weyl's criterion asserts that equidistribution is equivalent to having 
$$\sum_{a^2+1\equiv 0\pmod n}e\left(\frac{ak}{n}\right)=o\left(2^{\omega(n)}\right)$$
as $\omega(n)\rightarrow\infty$ for each $0<k<n$. By the chinese remainder theorem, the above statement equates to having
$$\prod_{p|n}\cos\left(\frac{2\pi k \overline{(n/p)}\rho_p}{p}\right)=o(1)$$
where $\rho^2_p+1\equiv 0\pmod p$ and $\overline{x}$ denotes the multiplicative inverse of $x \pmod p$. 
Since $0$ and $p/2$ cannot be roots $\pmod p$, the limits of such trigonometric products will  be zero if unboundedly many of the arguments do not converge to $0$ or $\pi$ as $\omega(n)\rightarrow\infty$, which amounts to having unboundedly many of the sequences 
$$\frac{n}{p}\pmod p\hspace{1cm}p|n$$
diverge. Is it possible that all but finitely many of them could converge? 
 A: Yes, this is a well-known. In fact one can prove the following: 
Let $\{\alpha_n\}$ be an arbitrary complex sequence with $|\alpha_n| \leq 1$ for all $n \geq 1$ and finite $\ell^2$-norm. Then for any positive numbers $D, N$ one has
$$\displaystyle \sum_{d \leq D} \sum_{v^2 + 1 \equiv 0 \pmod{d}} \left \lvert \sum_{n \leq N} \alpha_n e \left(\frac{vn}{d} \right)\right \rvert^2 \ll (D + N) \lVert \alpha \rVert_2^2,$$
where $\alpha = \{\alpha_n\}$ and $\lVert \cdot \rVert_2$ is the $\ell^2$-norm. 
The proof goes as follows. It clearly suffices to consider dyadic intervals of $D$, so we may assume that $D < d \leq 2D$ say. We then consider the equidistribution of the numbers $v/d$, with $v^2 + 1 \equiv 0 \pmod{d}, 0 \leq v < d$. For each such $v$ there exist positive integers $r,s$ such that $r^2 + s^2 = d$ and $r \equiv vs \pmod{d}$. This implies that we have an equation of the form
$$\displaystyle r - vs = \ell d,$$
which we rearrange as 
$$\displaystyle \frac{r}{sd} - \frac{\ell}{s} = \frac{v}{d}.$$
By exploiting symmetry we can always assume $|s| \geq |r|$, and that $s > 0$. Since $r^2 + s^2 = d$ this implies 
$$\displaystyle \frac{|r|}{sd} = \frac{|r|}{s(r^2 + s^2)} \leq \frac{|r|}{s(2|r|s)}   \leq \frac{1}{2s^2},$$
with equality only when $r = s = 1$. Hence if we have two pairs $(v_1, d_1), (v_2, d_2)$ say we get the equality
$$\displaystyle \frac{v_1}{d_1} - \frac{v_2}{d_2} = \frac{r_1}{s_1 d_1} - \frac{r_2}{s_2 d_2} - \frac{\ell_1}{s_1} + \frac{\ell_2}{s_2},$$
which implies that the mod 1 distance between $v_1/d_1, v_2/d_2$ is $\gg 1/(s_1 s_2) \gg 1/D$ (this is because $s_i \gg d_i^{1/2}$, and $d_1, d_2 \in [D,2D)$). The large sieve inequality then gives the conclusion above. 
With (a lot) more work one can even show that the fractions $v/p$ equi-distribute, where we take prime moduli only. See the following paper by Duke, Friedlander, and Iwaniec. 
