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