Using class field theory one can prove for every integer $n>0$ there exists a monic polynomial $f$ of degree $h(-4n)$ such that for any odd integer $m$ coprime to $n$ we have the following equivalence: 

$$\exists x,y\in\mathbb{Z}:(m=x^2+ny^2)\iff \exists x,y\in\mathbb{Z}:(f(x)\equiv 0\bmod m)\land (y^2\equiv-n\bmod m)$$

However counting the pairs $(x,y)\in\mathbb{Z}^2$ such that $m=x^2+ny^2$ can only be expressed in terms of similar "nice" sums of Jacobi symbols for just a few specific values of $n$. For the other cases you can simplify this to expressions that count the roots of polynomials modulo the primes dividing $m$. 

Though the formula you have for counting solutions to $m=x^2+y^2$ can be generalized, namely the number of proper representations of any integer $n$ by all the positive definite binary quadratic forms of some fixed discriminant can also be expressed as a similar divisor like sum involving Jacobi symbols using Dirichlet's Mass formula. While since the arbitrary representations of any $n\in\mathbb{N}$ by a quadratic form are proper representations of $n/d$ coordinate wise scaled by $\sqrt{d}$ for a square $d\mid n$. We can prove that if $r_2(n)=4(d_1(n)-d_3(n))$ is the number of integer pairs $(x,y)\in\mathbb{Z}^2$ which satisfy $x^2+y^2=n$ while if $s(n)=|\{d\mid n:\sqrt{d}\in\mathbb{N}\}|=\prod_{p\mid n}\left(1+\lfloor v_p(n)/2\rfloor\right)$ is the number of perfect squares dividing $n$ and $\omega(n)=\sum_{p\mid n}1$ is the number of distinct primes dividing $n$ then for every integer $m$ coprime to $n$ we define:

$$w(m)=\begin{cases}6&\text{ if }d=-3\\4&\text{ if }d=-4\\2&\text{ otherwise}\end{cases}$$
$$e_m(n)=\prod_{p\mid n}\frac{1}{2}\left(1+\left( \dfrac{m}{p}\right)\right)=\begin{cases}1&\text{ if }\exists x\in \mathbb{Z}:x^2\equiv m\bmod n\\0&\text{ if }\not\exists x\in \mathbb{Z}:x^2\equiv m\bmod n\end{cases}$$

So that if $f_1,\ldots ,f_{h(D)}$ are the reduced positive definite binary quadratic forms of discriminant $D$ where $h(d)$ is the class number of $D$ then for any odd integer $n\in\mathbb{N}$ coprime to $D$ we must have:

$$\small\sum_{k=1}^{h(D)}|\{(x,y)\in\mathbb{Z}^2:f_k(x,y)=n\land \gcd(x,y)=1\}|=w(D)\sum_{d\mid n}\left( \dfrac{D}{d}\right)=w(D)e_D(n)2^{\omega(n)}$$

$$\small\implies\sum_{k=1}^{h(D)}|\{(x,y)\in\mathbb{Z}^2:f_k(x,y)=n\}|=w(D)\sum_{d\mid n}\left( \dfrac{D}{d}\right)s(d)=w(D)e_D(n)\frac{r_2(n)}{4}$$

For instance if $D=-28$ then $f(x,y)=x^2+7y^2$ is the only reduced form of discriminant $D$ because $h(-28)=1$ therefore for every odd integer $n>0$ not divisible by $7$, when there exists integers $a,b\in\mathbb{Z}$ such that $n=a^2+7b^2$ then we must have by both of our previous formula that:


$$|\{(x,y)\in\mathbb{Z}^2:x^2+7y^2=n\land \gcd(x,y)=1\}|=2^{\omega(n)+1}$$
$$|\{(x,y)\in\mathbb{Z}^2:x^2+7y^2=n\}|=\frac{r_2(n)}{2}$$

For details on reduction theory read Chapter 4 of this: http://www2.math.ou.edu/~kmartin/ntii/ntii.pdf Whereas for Dirichlet's formula read page 5 of this: http://www2.math.ou.edu/~kmartin/ntii/chap4.pdf