Sum of squares diophantine equation

How can I count all integer solutions to the equation $$x^2+dy^2=n$$

given $n$ and $d$, where $n$ is composite (or prime)?

By counting I mean any algorithm faster than brute force enumeration.

What you want is the Hardy-Muskat-Williams algorithm, link HMW, which solves the Diophantine equation $a x^2 + c y^2 = n$ with the additional constraint that they give only the primitive representations, meaning $\gcd (x,y) = 1.$ If $n$ is not squarefree, you must solve the problem with the same $a,c$ but now $n/t^2$ for all such possible $t.$

It is no joke to program this. I have done it, largely because the first implementation in the language Mathematica was incorrect. It worked for $n$ prime but was unable to solve $x^2 + 14 y^2 = 65,$ which does have the solutions $x=\pm 3, \; y = \pm 2.$ The problem is that the article did not include any examples where some "seed" values (certain square roots $\pmod n$) failed to yield anything useful. So the Mathematica programmers assumed one seed value would be enough, whereas one must check all of them.

As you asked to count, just about the only example where you would be able to find an explicit formula that just counts the solutions, without finding them first, is the case $x^2 + y^2 = n.$ Such explicit counting formulas should be available for class number one, and maybe for the idoneal numbers, I am not sure in the latter case.

I assume $d$ is positive, else there may be infinitely many solutions.

I suspect there's nothing really fast. If you could find all the solutions to $x^2+y^2=n$ fast, then you could factor $n$ fast, and you'd be one wealthy mathematician. $x^2+dy^2=n$ has similar qualities.

EDIT: concerning factorization and counting solutions:

Suppose $r^2+ds^2=m$ and $t^2+du^2=n$. Think of $r^2+ds^2$ as the norm of $r+s\sqrt{-d}$ (or don't - it just makes it easier for me to see where the next formula comes from), then $(r+s\sqrt{-d})(t+u\sqrt{-d})=rt-dsu+(ru+st)\sqrt{-d}$ so $(rt-dsu)^2+d(ru+st)^2=mn$. So representations of a number are related to representations of its factors. This works out very simply in the case $x^2+y^2$, where you can derive a formula for the number of representations of $n$ in terms of the prime divisors that are congruent to 1 modulo 4. It can be more messy for other values of $d$ and I don't know it well enough to write it out on the spot and I don't have access to my references. But there are plenty of books and websites that talk about representations by positive definite binary quadratic forms, so the information is out there, you just have to hunt for it.