Efficient software implementation of $x^2+3y^2=N$ I would like to implement a solver (in C) for the Diophantine equation $x^2+3y^2=N$ for non-negative integers $\{x,y\}$ and positive integer $N$. I have read online that one has to prime factorize N and express the primes in $\mathbb{Z}[\sqrt{-3}]$. Unfortunately I do not have sufficent mathematical expertise to follow through the heavy proofs for the last part.
Is there an easy-to-understand method of expressing primes in $\mathbb{Z}[\sqrt{-3}]$ that can also be efficiently implemented in software?
PS: My background in mathematics consists of engineering calculus and basic number theory up to modular arithmetic.
 A: I guess in your case it might be useful to learn about Eisenstein integers. They are given by
$$ z= a+ b \omega$$ with $a$ and $b$ integers and $\omega = \exp(2\pi i/3) = (i\sqrt{3} -1)/2$.
Using the Eistenstein integers, your equation can be factorized
$$x^2 +3 y^2 = \Bigl(x+ (1+2 \omega) y  \Bigr)\Bigl(x+ (1+2 \omega^2)y\Bigr)= N.$$
So given the unique prime-factorization of Eisenstein numbers (up to the units $\pm 1,\pm \omega,\pm \omega^2$), we have that $N$ should be factorizable in the same form.
So you can first factoring $N$ in the ordinary integers and thus we only need to discuss the case where $N$ is a prime number $p$. Now there are two cases
1) $p = -1 \pmod 3$: then $p$ is an Eisenstein prime and the equation has no solution (this is the case if any of the prime factors of $N$ are of this form).
2) in the other cases, we can factoring $p = (a + b \omega)(a+ b \omega^2)
= a^2-ab +b^2$, we can then compare the prime factors to determine $x$ and $y$.
So the problem for the numerical algorithm is to solve $p=a^2-ab+b^2$ for $p$ a prime congruent to $1$ modulo 3 for which you know that a solution exists (the cases $p=3$ is covered by $a=1$, $b=2$). I have to admit that I do not know about any efficient algorithm for that except for trying $a$ and $b$.
Edit:
here, is an algorithm which given a prime $p$ which is congruent to 1 modulo 3 finds $a$ and $b$ with $p= a^2 -a b +b^2$.
First, we note that $\gcd(a,b)=1$, that is $a$ and $b$ are coprime. All coprime integers can be generated by the Farey sequence where Wikipedia has an imlementation in Python. If the prime $p$ is not too large the algorithm of taking the coprime integers from the Farey sequence and simply testing if $p=a^2-ab+b^2$ will be efficient enough.
A: Although the equation is simple. But it often occurs in many cases. Required to represent the sum of the squares of the product of the multipliers. And the inverse problem. To represent the number as a sum.
Therefore, it is convenient equation to represent and solve it like this:
$$x^2+3y^2=N=ab$$
Then the solution can be represented as:
$$x=(p-s)kn+(p^2+5ps)n^2$$
$$y=pkn+(p^2-2ps+2s^2)n^2$$
$$a=k^2+2pkn+(p^2+12s^2)n^2$$
$$b=(4p^2-2ps+s^2)n^2$$
Or so:
$$x=k^2-2(4p+s)kn+(7p^2+26ps-8s^2)n^2$$
$$y=k^2-4(2p-s)kn+(7p^2-16ps+4s^2)n^2$$
$$a=4(k^2-(2p-s)kn+(p^2-ps+7s^2)n^2)$$
$$b=(k-(7p-2s)n)^2$$
$p,s,k,n$ - integers of any sign.
It is necessary to consider not mutually simple solutions. Or Vice versa. Reduce to obtain relatively simple solutions.
