Algorithm to find solutions $(p,x,y)$ for the equation $p=x^2 + ny^2$.

Question

As the classical book of David Cox argues,

Assume the conditions are satisfied and $p$ can be represented as $x^2 + ny^2$. What would be a way to find solutions $(p,x,y)$ efficiently? Ideally, one would have one algorithm that works for all $n$.

Answer 1

solve $$ \beta^2 \equiv -4n \pmod p. $$ by Cipolla or Tonelli-Shanks. If $\beta $ is odd, replace by $p-\beta,$ so now we have an even number, call it $b,$ with $$ b^2 \equiv -4n \pmod {4p}. $$ That is $$ b^2 = -4n + 4 p t,$$ and $$ b^2 - 4 p t = -4n. $$ So we have a positive binary quadratic form $$ \langle p, b, t \rangle $$ of discriminant $-4n.$ Find the Gauss reduction, you give the hypothesis that this is guaranteed to be $$ \langle 1, 0, n \rangle. $$ Gauss reduction is very fast, a step by step process that is comparable to the Euclidean algorithm for gcd, both in speed and style of the steps.

The inverse of the reduction matrix shows how to get from $ \langle 1, 0, n \rangle $ back to $ \langle p, b, t \rangle , $ including the representation of $p.$

Give names

$$ P = \left( \begin{array}{rr} p & \frac{b}{2} \\ \frac{b}{2} & t \end{array} \right) $$

and

$$ N = \left( \begin{array}{rr} 1 & 0 \\ 0 & n \end{array} \right). $$

The thing I call a 'reduction matrix' above is a 2 by 2 matrix $R$ of integers of determinant $+1,$ with $$ R^T P R = N. $$

Let $$ Q = R^{-1}, $$ this is easy to find for the 2 by 2 case with determinant 1. Then $$ Q^T N Q = P. $$ If you write it out with explicit entries, you will see that the left column of $Q$ gives the representation of $p$ as $x^2 + n y^2.$

I recommend you go through this process to find out how the primes $$ 59, 101, 167, 173, 211, 223, 271, \ldots $$ are represented by $x^2 + 23 y^2.$ in this case, the polynomial in Theorem 5.1 above is $$ f(z) = z^3 - z + 1. $$