PARI code analyzing whether given integer is presentable as sum of two integer squares

[victor miller]

Alex, Here's a post from Noam Elkies which gives some pari code of
mine.  This works if p = 1 mod 4.  If you have a composite you need to
factor it (at least into prime powers for primes = 1 mod 4) to get all
solutions.

Victor


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

In article , rlr 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.



On Nov 21, 2007 8:03 PM, Alexander Povolotsky <apovolot at gmail.com> wrote:
> Hi,
>
> Could anyone write/post PARI code (function) analyzing whether given
> integer is presentable as sum of two integer squares (finding all the
> combinations if more than one exist)?
> I didn't come across of such code currently posted in any OEIS related
> sequence ...
>
> Thanks,
> Alex
>