What this does is to use the well known construction (I think that it's in Hardy and Wright), that says that if $0 < a < (p-1)/2$ satisfies $a^2 = -1 \bmod{p}$ if you run the continued fraction algorithm for $a/p$ "half-way" to get the convergent $r/s$ then $r^2 + s^2 = p$. What the above one-liner does is to set up the lattice $$\pmatrix {a&p \cr p 1 & 1\cr}$$ 0 \cr}$$The shortest vector in this lattice has L^2 norm of p. Post Undeleted by Victor Miller 2 improved formatting Thanks to Noam Elkies for telling me to post my one-liner to solve this in gp (this dates from around 1993): fermat(p) = qflll([lift(sqrt(Mod(-1,p))),p;1,0])[1,] What this does is to use the well known construction (I think that it's in Hardy and Wright), that says that if 0 < a < (p-1)/2 satisfies a^2 = -1 \bmod{p} if you run the continued fraction algorithm for a/p "half-way" to get the convergent r/s then r^2 + s^2 = p. What the above one-liner does is to set up the lattice$$\left(\begin{matrix} $\pmatrix {a&p \cr p & 1 \end{matrix}\right)$$[for some reason the matrix environment isn't working correctly]. 1\cr}$$ The shortest vector in this lattice has$L^2$norm of$p\$.