Expressing primes $p\equiv 1 \pmod 3$ in the form $p = x^2 + xy + y^2$ 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 explain that these are the primes $\not\equiv 2\mod 3$), as well as an efficient algorithm to obtain such a factorization given $p$.

*

*Shoup, Victor, A computational introduction to number theory and algebra, Cambridge: Cambridge University Press (ISBN 978-0-521-51644-0/hbk). xvii, 580 p. (2009). ZBL1196.11002.

 A: 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 or Schoof's algorithm.
