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

S(D,x) = { q < x : D is a quadratic residue mod 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!