Let D$D$ be a fixed positive squarefree integer. For a positive integer x$x$, define
S(D,x) = { q < x : D is a quadratic residue mod q }$S(D,x) = \{ q < x : D \text{is a quadratic residue} \pmod q \}$.
Here q$q$ can be any integer, not necessarily a prime. Are elements of S(D,x)$S(D,x)$ evenly distributed? In other words, let 0 < a < b < 1$0 < a < b < 1$ be constants and consider an interval I = [ax,bx]$I = [ax,bx]$. Ideally I would like to see some result that says that the number of elements of S(D,x)$S(D,x)$ in I$I$ is proportional to the length of I$I$ on the average as x$x$ goes to infinity (is this true?).
I am aware of classical 1917-1918 results of Vinogradov and Polya (and some later developments) about distribution of quadratic residues, which in particular imply that quadratic residues modulo a fixed prime p$p$ are evenly distributed in the interval [0,p]$[0,p]$ in the same sense as I described above. What I need, however, is a result on the distribution of moduli with respect to which a fixed integer is a quadratic residue, and I cannot find anything like this in the literature.
In other words, I am wondering how the divisors of x^2-D$x^2-D$ are distributed on the average as x$x$ goes to infinity. It is a well-known fact that for an arbitrary integer y$y$, there are unproportionally many (on the average) small and large divisors of y$y$ as y$y$ goes to infinity, i.e., divisors are not uniformly distributed. But what if we take y$y$ in the special form x^2-D$x^2-D$?
Any thoughts on the subject are much appreciated!