Find efficient way to generate all solutions to Diophantine equation $a^2+5ab+3b^2-c^2=0$ under a given bound $N$ I am looking to solve Diophantine equation $a^2+5ab+3b^2-c^2=0$.
a, b, c are all positive.
Since the number of solutions are infinite. Lets say we are only interested in solutions till a limit N ie $1 \le a, b, c \le N$. If you have a expression for solution, I would want it to be better than trying all values of a, b till the limit $N$.
We need all solutions set (not the number of solutions).
I already have a solution which reads like this:

Consider rational solutions to:
$$x^2+5xy+3y^2=1.\tag1$$
We have $(-1,0)$ is a rational solution.
If $p,q$ are relatively prime integers, there are usually two values
$t$ such that $(-1+pt,qt)=(-1,0)+t(p,q)$ is a solution, one with
$t=0.$
$$\begin{align}0&=(-1+pt)^2+5(-1+pt)qt+3(qt)^2-1\\&=-2pt+p^2t^2-5qt+5pqt^2+3q^2t^2\\&=t(\left(p^2+5pq+3q^2)t-(2p+5q)\right)
> \end{align}$$
So $t=0$ or $t=\frac{2p+5q}{p^2+5pq+3q^2}.$
The.
$$(x,y)=\left(\frac{p^2-3q^2}{p^2+5pq+3q^2},\frac{2pq+5q^2}{p^2+5pq+3q^2}\right)$$
And you get $$(a,b,c)=(p^2-3q^2,2pq+5q^2,p^2+5pq+3q^2).\tag2$$
If $(p,q)=(2,1)$ you get your solution $(1,9,17).$
This won't give primitive solutions, in general. When $(p,q)=(5,-2)$
then $(a,b,c)=(-13,0,-13),$ for example.
If $\gcd(p,q)=1,$ we can show $\gcd(a,b)=1$ or $13.$ It can only be
$13$ if $p\equiv 4q\pmod{13}.$

Here we have reduced the equation in $a,b,c$ to $p,q$. For any value of p,q we would get a solution. However it is still not very computationally efficient. Lets say we have to generate all such solutions till $N$ ie $a,b,c \le N$ .
If I iterate over all values of $-N < p \le N$ and $-N < q \le N$ which makes it's net computational cost $O(N^2)$. It still doesnt guarantee to find all solutions within N. It may miss a few and also it may generate a few false positive (ie >N) . We don't really know which are those p,q which would yeild all solutions which would be $(a,b,c)<= N$.  Can you suggest improvements to the solution. Ideally this should be computed in time proportional to the number of solutions and we need not do any checking.
PS: this is a followup question to Solve Diophantine equation $a^2+5ab+3b^2-c^2=0$
The approached for this versions of the questions where are interested in bounded solution are going to be different from the original one and hence should not be merged with the original question.
Assumption of $-N \le p,q \le N$ is over-simplification
Obviously we are looking for more tighter bound of $p,q$ if we solve

*

*$1 \le p^2-3q^2 \le N$


*$1 \le 2pq+5q^2 \le N$


*$1 \le p^2+5pq+3q^2 \le N$
Which is what this question is all about.
If you have any suggestion to take it to anything better than $O(N^2)$ you are welcome.
You can assume we have $O(N^2)$ memory and arithmetic operations +-*/% as $O(1)$
 A: Since this can be written as:
$$(2a+5b)^2-13b^2=(2c)^2,$$ we will first interest ourselves in solutions to $u^2-13v^2=w^2$ which can be written as $$(u-w)(u+w)=13v^2.$$
Ultimately, we'll want $u\equiv v\pmod 2$ and $w$ even.
Writing $v=v_0v_1v_2$ where $v_0$ is square-free, and $$u+w=13v_0v_1^2, u-w=v_0v_2^2,$$ we get:
$$u=\frac{v_0(13v_1^2+v_2^2)}{2}\\
w=\frac{v_0(13v_1^2-v_2^2)}2$$
Now, $w$ can only be even if $v_1,v_2$ are both even, or both odd. (Since $v_0$ is square-free, it can't be divisible by $4.$)
Also, $u\equiv v\pmod 2$ only if $v_1,v_2$ are both even, or if $v_1,v_2$ are both odd.
Now, if $(a,b,c)$ is a primitive root, then $\gcd(a,b)=1.$ But then $\gcd(2a+5b,b)=\gcd(2,b),$ so we can restrict to the case when $v_0=1,2,$ and $\gcd(v_1,v_2)=1,2.$ because any other factor of $v_0$ gives a common factor.
Then $$\begin{align} 
b&=v_0v_1v_2, \\
u&=\pm\frac{v_0(13v_1^2+v_2^2)}{2}\\
a&=\frac{u-5b}2=\frac{v_0(\pm(13v_1^2+v_2^2)-10v_1v_2)}{4},\\
c&=\left|\frac{v_0(13v_1^2-v_2^2)}4\right|
\end{align} $$
Where the $\pm$ for $a$ is the same as the $\pm$ for $u.$
Note we again get common factors of $13$ when $13\mid v_2.$ But we get the same tuple divided by $13$ if we compute for $v_1'=v_2/13, v_2'=v_1.$
But if $v_0=2$ we always get a common factor between $a,b.$
So, we iterate through all pairs $v_1,v_2$ of the same parity with $\gcd(v_1,v_2)=1,2$ and $|v_1v_2|\leq N$ and $v_2$ not divisible by $13.$ There are $O(N\log N)$ such pairs. If we cache the prior values for the GCD, we can iterate over these in $O(N\log N)$ time.
Some of $a,c$ will be outside the range. We can skip these.
Note, if $13\mid v_2,$ the primitive solution we get related to the above approach is by dividing by $13,$ and thus the values become smaller. We can handle that with some trickery. Essentially, if $13\not\mid v_1,$ you can solve for $v_1,v_2$ and for $v_1'=v_1,v_2'=13v_2.$
If we want to generate non-primitive solutions, we take our primitive solution by multiplying by $k$ for $k\leq \frac{N}{\max(|a|,|b|,|c|)}.$
I'll try to write a Python script when I get to my computer.
