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 involving root counts 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$ whereas 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 if for every integer $m$ coprime to $n$ we define:
$$w(m)=\begin{cases}6&\text{ if }m=-3\\4&\text{ if }m=-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}$$
We see if $\mathscr{F}_D$ is the set of reduced positive definite binary quadratic forms of discriminant $D$ so that by definition $h(D)=|\mathscr{F}_D|$ then for an arbitrary odd integer $n\in\mathbb{N}$ coprime to $D$ we must have that:
$$\small\sum_{f\in \mathscr{F}_D}|\{(x,y)\in\mathbb{Z}^2:f(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_{f\in \mathscr{F}_D}|\{(x,y)\in\mathbb{Z}^2:f(x,y)=n\}|=w(D)\sum_{d\mid n}\left( \dfrac{D}{d}\right)s(d)=w(D)e_D(n)(d_1(n)-d_3(n))$$
Thus in particular when $D=-4$ since $f(x,y)=x^2+y^2$ is the only reduced form of discriminant $D$ these formula simplify to your original identity for counting the number of representations of an integer as a sum of two squares. While note 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\}|=2(d_1(n)-d_3(n))$$
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