From: Noam D. Elkies <elkies@mafiveth.harvard.edu>
Subject: Re: Decomposing a Legendre Prime quickly?
Date: 16 Apr 2000 13:51:24 GMT
Newsgroups: sci.math.research
Summary: [missing]

In article <8dashl$2cp$1@news1.abac.com>, rlr  <rlr@earth.sun.galaxy> wrote:

>Can anyone give an efficient procedure for finding the 2 integers whose
>square sum is a Legendre prime? A Legendre prime is a prime of form 1 mod 4.

\begin{rant}[mild,terminology]
"Legendre prime"?  I don't think this is standard terminology; the theorem
that an odd prime p can be written as a^2+b^2 if and only if p is 1 mod 4,
in which case the representation is unique, is usually attributed to Fermat,
but "Fermat prime" has a different meaning, so something like "4m+1 prime"
is probably best -- it's unambiguous, and "4m+1" is shorter than either
"Fermat" or "Legendre"...
\end{rant}

\begin{answer}
Anyway, to answer your question: in late 1992, Victor Miller posted
the following one-liner:

  fermat(p) = lllint([lift(sqrt(mod(-1,p))),p;1,0])[1,]

this is a GP program; having entered this into GP one can write

  fermat(200420221)

and get

  [-10011, -10010]

in reply.  The program uses several sophisticated built-in GP functions, but
the basic idea is: find x such that p|x^2+1 (this is the "sqrt(mod(-1,p))"
in Victor's program), and expand x/p in a continued fraction to find minimal
a,b such that a/b is x mod p; then p=a^2+b^2.  To find x, choose nonzero
c mod p at random, and compute c^m mod p (recall that p=4m+1); half the time
(if c is a quadratic nonresidue of p), this will be a square root of -1 mod p.
\end{answer}

\begin{sig}[spamblock_clue]
--Noam D. Elkies (remove Fermat prime from e-address to reply)
  Department of Mathematics, Harvard University
\end{sig}