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$.