Statement of the problem
Firstly, consider the infamous Pell's equation:
$x^{2}-dy^{2}=1$. Here $x$ and $y$ are integers and $d$ is a nonsquare integer. It is known ([3]) that all solutions of this equation can be obtained from the solution with least positive $x$ called the fundamental solution. By the way, it is convenient to write solutions in the form $a+\sqrt{d}b$. This way if we denote the fundamental solution by $x_{1}+\sqrt{d}y_{1}$ then the $k$-th solution will be $x_{k}+\sqrt{d}y_{k} = (x_{1}+\sqrt{d}y_{1})^{k},k\geq 1$. And I need to note here that this equation is always solvable because we can obtain the fundamental solution from the continued fraction expansion of $\sqrt{d}$.
Now, let's look at the generalised Pell's equation: $x^{2}-dy^{2}=N$. Here $x$,$y$ and $N$ are integers and $d$ is a nonsquare integer. It turns out ([1],[2]) that the solutions (if it's solvable, cause it may not have a solution) all have quite a simple structure:
$
\eta_{m}\epsilon_{0}^{k}, 1\leq m\leq h,k \geq 0.
$
Where $\eta_{m}$ is the fundamental (with least positive $x$ component) solution of $m$-th equivalence class and $\epsilon_{0}$ is the fundamental solution of standard Pell's equation $x^{2}-dy^{2}=1$.
I'm talking about equivalence classes so here is the definition of equivalence: two solutions $f_{1} = a_{1}+\sqrt{d}b_{1},f_{2} = a_{2}+\sqrt{d}b_{2}$ of the generalised Pell's equation are equivalent if their ratio is a solution of the standard Pell's equation. Note that it doesn't matter if we take $\frac{f_{1}}{f_{2}}$ or $\frac{f_{2}}{f_{1}}$. Because solving the Pell's equation is essentially finding numbers of required norm in the ring $\mathbb{Z}+\sqrt{d}\mathbb{Z}$.
Finally, let's denote by $f(N,d)$ that number $h$ of equivalence classes of solutions to the generalised Pell's equation. This is what I want to research. Experimental results are as follows (changing $N$,$d$ is fixed): it behaves quite randomly, but the average seems to tend to a constant. This motivates me to consider the function: $$ F(M,d) = \sum_{N = 1}^{M} f(N,d) $$ My aim here is to get a good approximate formula for this function and then take a limit $\lim_{M \rightarrow +\infty} \frac{F(M,d)}{M}$ to see if it really tends to a constant.
My current approach
If we could somehow get the Dirichlet series $g(s) = \sum_{N=1}^{+\infty} \frac{f(N,d)}{N^{s}}$ then we'd just use Perron's inversion formula to recover the sum. The hardest thing for me is that $f(N,d)$ doesn't seem to have any simple properties that I could use to get some functional equations, etc. It feels like pure randomness and I'm stuck looking from this angle. For instance, trying $\sum_{x,y >0} \frac{1}{(x^{2}-d y^{2})^{s}(x+\sqrt{d}y)^{s-1}}$ yields something that, roughly speaking, behaves like a simplified version of that Dirichlet series near $s = 1$ but summing it is problematic for me. Because we have a double sum and this means the error terms should be quite small to be able to sum it properly.
My questions
- How difficult is this problem?
- Are there any tools known (maybe from analytic number theory) to work with such problems?
Some literature
- Matthews K.R.: The diophantine equation $x^2-Dy^2=N, D>0.$ Expositiones Mathematicae, 18, 323-331 (2000)
- Mollin R.A.: Simple Continued Fraction Solutions for Diophantine Equations. Expositiones Mathematicae, 19, 55-73 (2001)
- Clark, Pete. Number Theory: A Contemporary Introduction. University of Georgia. 93-102
