Integer as sum of two squares of integers

[tytus]

How to efficently check if given integer (from range 1-10^12) can be represented as sum of two squares of integers, so if n=a^2+b^2.
I know that:
a) Those numbers n whose prime divisors of the form 4*m+3 divide n an even number of times can be written as the sum of two squares.
b) Jacobi's Two Square Theorem: The number of representations of a positive integer as the sum of two squares is equal to four times the difference of the numbers of divisors congruent to 1 and 3 modulo 4.

But when n is big (for example 10^10 or 10^12) it's hard to find all prime divisors of form 4*m+3 or all divisors congruent to 1 and 3 modulo 4.

Can anyone give me a hint how to efficently check if n can or cannot be represent as sum of two squares of integers??

[wblipp]

Originally Posted by tytus

But when n is big (for example 10^10 or 10^12) it's hard to find all prime divisors

http://www.alpertron.com.ar/ECM.HTM and dozens of other programs will factor such numbers in well under a second.

[tytus]

OK - I see... But I would like to do it myself... I mean - if I would like to write my own solution, so how can I check efficently if number n can be represented as sum of two squares of integers?? Can anyone give me a hint??

[ewmayer]

Originally Posted by tytus

OK - I see... But I would like to do it myself... I mean - if I would like to write my own solution, so how can I check efficently if number n can be represented as sum of two squares of integers?? Can anyone give me a hint??

As a simple start, to see if n = x^2 + y^2, just loop over candidate integers x (presieved along the lines you mention), calculate z = n - x^2 (for n = O(10^12) you'll want to cast x to double or 64-bit int before squaring and use the same type for z), and for each x check if z is a perfect square via e.g. (using C notation - and I haven't tested this, so no warranty ;))

int x;
double n, y, z;
/* floating-point fudge factor; assumes sqrt()
   implementation has RO error <= this value */
const double eps = 1.0e-12;
...
{init n}
...
(loop over x, assume x increases with each loop execution)
{
    z = (1.0*x);
    z = n - z*z;
    if(z < 0)
        break;
    /* assume FP sqrt inexact, so add small fudge factor before truncating */
    y = floor(sqrt(z) + eps);
    if(z == y*y)
        printf("n = %d^2 + %d^2\n", x, (int)y);
}

For n = O(10^12) one only needs to test x up to around a million, which should take almost no time at all even on a modest PC. If speed is an issue there are added things one can use to reduce the number of such tests, e.g. a fast floor-less NINT that uses the well-known add-and-then-subtract-floating-rounding-constant trick, and so forth.

Have fun,
-Ernst

[philmoore]

There is an algorithm for expressing a prime p = 1 mod 4 as the sum of two squares, which I have used in the past, but need to review before I can post details. I just remember that it starts by finding a quadratic non-residue. I suspect that Dario Alpern ("alpertron" on this forum) uses the same routine in his Java applet. I'll pm him for details, otherwise, I'll eventually try to get the details posted here myself.

[alpertron]

Yes, I'm using that routine. See this page for details.

This means that for composite numbers, we need to factor it. Notice that we don't need to factor the number if we want to express it a sum of three or four squares. See my Sum of squares applet where you can enter expressions. You can even see the decomposition on squares of RSA-2048 with this applet.

[philmoore]

Thanks, Dario, for the link to your web-page. You are a teacher as well as a scholar! Your algorithm looks like Fermat's method of descent, which is also the first method I used to solve this problem. I eventually found another algorithm called "Cornacchia's algorithm" which solves the problem with fewer steps. I haven't had a chance to look at my notes until tonight or I would have posted earlier. (Having a 5-month old baby certainly seems to limit my free time!) So here is the algorithm:

p is a prime = 1 mod 4
Solve p = x^2 + y^2:

First find a quadratic non-residue b to p. (As Dario has said, the easiest way is to check all prime numbers in sequence until we find b with Jacobi(b,p)=-1.)
x = b^((p-1)/4) mod p
sqrtp = floor(sqrt(p))
y = p%x
if y=1 then done else
do while x>sqrtp
{r = x%y
x = y
y = r}

[Raman]

Originally Posted by philmoore

p is a prime = 1 mod 4
Solve p = x^2 + y^2:

First find a quadratic non-residue b to p. (As Dario has said, the easiest way is to check all prime numbers in sequence until we find b with Jacobi(b,p)=-1.)
x = b^((p-1)/4) mod p
sqrtp = floor(sqrt(p))
y = p%x
if y=1 then done else
do while x>sqrtp
{r = x%y
x = y
y = r}

Wondering why this algorithm does work out?

By the way, for the extension, what is the modification needing to be done
to write out the given prime
p = 1, 3 (mod 8) into x²+2y²
p = 1 (mod 3) into x²+3y²

p = 1, 9 (mod 20) into x²+5y²
forms?

[LaurV]

Originally Posted by Raman

..... write out the given prime p = 3 (mod 8) into x²+2y²

Is this always possible? Is there any efficient algorithm doing this? Can you point me to some theory paper related to this? Thanks in advance.

[Raman]

Originally Posted by LaurV

Is this always possible? Is there any efficient algorithm doing this? Can you point me to some theory paper related to this? Thanks in advance.

I checked through to the values for p = 1000, the answer is being always yes!

p = 1 (mod 4) uniquely writable into x²+y²
p = 1, 3 (mod 8) uniquely writable into x²+2y²
p = 1 (mod 3) uniquely writable into x²+3y²
what else?
p = 1, 9 (mod 20) uniquely writable into x²+5y²
...
the forms!
cleanly, clearly

http://en.wikipedia.org/wiki/Fermat%...of_two_squares

(a²+Nb²)(c²+Nd²) = (ac+bd)²+N(ad-bc)² = (ac-bd)²+N(ad+bc)²

The algorithm please?
I think we need to find out the two squares (mod p),
such that x²+ky² = 0 (mod p), and then reduce it very similar
to the GCD process, so, thus

[Gammatester]

Originally Posted by LaurV

Is this always possible? Is there any efficient algorithm doing this? Can you point me to some theory paper related to this? Thanks in advance.

Use Cornacchia's algorithm (see http://en.wikipedia.org/wiki/Cornacchia%27s_algorithm or Crandall/Pomerance, Algorithm 2.3.12), it is very effective.