Let $D$ be a fixed positive squarefree integer. For a positive integer $x$, define

$S(D,x) = \{ q < x : D \text{is a quadratic residue} \pmod q \}$.

Here $q$ can be any integer, not necessarily a prime. Are elements of $S(D,x)$ evenly distributed? In other words, let $0 < a < b < 1$ be constants and consider an interval $I = [ax,bx]$. Ideally I would like to see some result that says that the number of elements of $S(D,x)$ in $I$ is proportional to the length of $I$ on the average as $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$ are evenly distributed in the interval $[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$ are distributed on the average as $x$ goes to infinity. It is a well-known fact that for an arbitrary integer $y$, there are unproportionally many (on the average) small and large divisors of $y$ as $y$ goes to infinity, i.e., divisors are not uniformly distributed. But what if we take $y$ in the special form $x^2-D$?

Any thoughts on the subject are much appreciated!