The Math Forum

Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Math Forum » Discussions » sci.math.* » sci.math

Topic: Fermat's theorem on sums of two squares (computational complexity

Replies: 8   Last Post: Aug 16, 2011 7:24 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Bill Dubuque

Posts: 1,739
Registered: 12/6/04
Re: Fermat's theorem on sums of two squares (computational complexity question)
Posted: Aug 16, 2011 2:07 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

Phil Carmody <> wrote:
> Pubkeybreaker <> writes:
>> On Aug 11, 7:27Ž?am, David Bernier <> wrote:
>>> Concerning Fermat's sum of two squares theorem, that every
>>> prime p of the form p = 4n+1 can be expressed as
>>> the sum of two perfect squares, books and on-line
>>> references don't say much about the computational
>>> complexity of decomposing a generic prime of the form
>>> 4n+1 as a sum of two squares.
>>> The PARI-gp package finds that the 51-digit number N51:
>>> N51 = 314,159,265,358,979,323,846,264,338,327,950,288,419,716,939,937,673
>>> is a prime p of the form 4n+1
>>> Is it feasible today to have a decent chance
>>> of decomposing a prime such as N51 as a sum of two squares ?

>> Yes. It is trivial. Look up "Cornacchia's Algorithm"

> Sometimes Cornacchia-Smith. It seems as if Smith may have
> had the idea earlier, but they are as near as dammit.
> In pari, the cheat or cheap method to decompose
> such a prime P into 2 squares is to factor P+0*I.
> The answer will come back (a+b*I)*(a-b*I), so a^2+b^2=P.
> I have no idea what algorithm it uses, so it may not
> be efficient asymptotically, but I've never seen it
> break a sweat (OK, it does a probable primality test
> first, but that's not so expensive).

This can be done simply by using the Euclidean algoritm in Z[i] to
compute gcd(p,i-j) = a+bi, where j^2 = -1 (mod p). Then p = a^2+b^2.
Below is a Macsyma program that does such. See also my posts here

/* For p=1 (mod 4) compute w = a+bi with p = ww' = a^2 + b^2 */

psplit(p) := block([j:0, e:floor(p/4), algebraic:true, x, y ],

block([modulus:p], /* search for j = sqrt(-1) (mod p) */
for a:2 while j^2 # -1
do j:rat(a)^e ),

x:p, y : j - %i, /* find gcd(p,j-i) in Z[i] by Euclidean algorithm */

do ( if y = 0 then return(x),
j : rat(x - y * floor(rat(x/y))),
x:y, y:j ))$

--Bill Dubuque

Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© The Math Forum at NCTM 1994-2017. All Rights Reserved.