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 $$\pmatrix {a&p \cr p & 1\cr}$$ The shortest vector in this lattice has $L^2$ norm of $p$.

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 $$\pmatrix {a&p \cr 1 & 0 \cr}$$ The shortest vector in this lattice has $L^2$ norm of $p$.

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} a&p \\ p & 1 \end{matrix}\right)$$ [for some reason the matrix environment isn't working correctly]. The shortest vector in this lattice has $L^2$ norm of $p$.

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 $$\pmatrix {a&p \cr p & 1\cr}$$ The shortest vector in this lattice has $L^2$ norm of $p$.