Fermat famously showed that the only primes $p$ of the form $x^2 + y^2$ are the primes such that $p \equiv 1 \mod{4}$. Furthermore, we now know “effective” versions of Fermat's theorem, i.e. given a prime $p$ such that $p \equiv 1 \mod{4}$, we know how to find integers $x$, $y$ such that $x^2 + y^2 = p$ in time polynomial in $\log p$ (see, for example section 4.5 in [1]). I would like an analogous theory for primes of the form $x^2 + xy + y^2$. In other words, I would like a precise characterization of which primes $p$ can be expressed in this form (**EDIT:** The [comments](https://mathoverflow.net/questions/371622/expressing-primes-p-equiv-1-pmod-3-in-the-form-p-x2-xy-y2#comment939751_371622) explain that these are the primes $\not\equiv 2\mod 3$), as well as an efficient algorithm to obtain such a factorization given $p$.


1. *Shoup, Victor*, [A computational introduction to number theory and algebra](https://www.cambridge.org/core/books/computational-introduction-to-number-theory-and-algebra/8C43512C31D9FDE058E360441175966C), Cambridge: Cambridge University Press (ISBN 978-0-521-51644-0/hbk). xvii, 580 p. (2009). [ZBL1196.11002](https://zbmath.org/?q=an:1196.11002).



Expressing primes $p\equiv 1 \pmod 3$ in the form $p = x^2 + xy + y^2$
This is an elaboration of the answer that Noam Elkies provided in the comments.


Suppose that $p=x^2 + xy + y^2$. Then note that $x$ and $y$ are small relative to $p$ (at most half as many digits). Note also that if $\zeta \not\equiv 1\pmod p$ satisfies $\zeta^3 \equiv 1\pmod p$ then $\zeta^2 + \zeta + 1 \equiv 0 \pmod p$, so
$$(x - \zeta y)(x - \zeta^2 y) = x^2 - (\zeta+\zeta^2)xy + \zeta^3 y^2 \equiv x^2 + xy + y^2 \equiv p \equiv 0 \pmod p.$$
Therefore either $x \equiv \zeta y \pmod p$ or $x \equiv \zeta^2 y \pmod p$; in the latter case we have $\zeta x \equiv y \pmod p$. This means that in the 2-dimensional integer lattice generated by the vectors $(1,\zeta)$ and $(0,p)$, there is an unusually short vector $(y,x)$ or $(x,y)$, which can be found by lattice-basis reduction as long as we have $\zeta$.


It remains to find $\zeta$. Formally, we can write
$$\zeta := {\sqrt{-3} - 1 \over 2},$$
and it is easy to check that if we can find a square root of $-3$ modulo $p$ then this formula does indeed give us a cube root of unity modulo $p$. But computing the square root can be done using the [Tonelli–Shanks algorithm](https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm) or [Schoof's algorithm](https://math.stackexchange.com/q/1999286/502981).


