Search All of the Math Forum:

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

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

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

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

 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

Phil Carmody <thefatphil_demunged@yahoo.co.uk> wrote:
> Pubkeybreaker <pubkeybreaker@aol.com> writes:
>> On Aug 11, 7:27?am, David Bernier <david...@videotron.ca> 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
http://math.stackexchange.com/questions/5877

/* 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

Date Subject Author
8/11/11 David Bernier
8/11/11 Pubkeybreaker
8/11/11 David Bernier
8/13/11 David Bernier
8/13/11 gnasher729
8/14/11 Phil Carmody
8/16/11 Bill Dubuque
8/16/11 Phil Carmody
8/16/11 David Bernier