mersenneforum.org Integer as sum of two squares of integers
 Register FAQ Search Today's Posts Mark Forums Read

 2005-02-12, 14:56 #1 tytus   Jan 2005 416 Posts Integer as sum of two squares of integers 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??
2005-02-12, 17:17   #2
wblipp

"William"
May 2003
Near Grandkid

2·1,187 Posts

Quote:
 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.

 2005-02-12, 18:23 #3 tytus   Jan 2005 410 Posts 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??
2005-02-22, 05:06   #4
ewmayer
2ω=0

Sep 2002
República de California

1175510 Posts

Quote:
 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 ;))

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

Last fiddled with by ewmayer on 2005-02-22 at 18:02

 2005-02-22, 19:24 #5 philmoore     "Phil" Sep 2002 Tracktown, U.S.A. 25×5×7 Posts 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.
 2005-02-22, 19:45 #6 alpertron     Aug 2002 Buenos Aires, Argentina 101110110012 Posts 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.
 2005-02-24, 07:22 #7 philmoore     "Phil" Sep 2002 Tracktown, U.S.A. 25×5×7 Posts 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}
2012-03-14, 17:35   #8
Raman
Noodles

"Mr. Tuch"
Dec 2007
Chennai, India

3·419 Posts

Quote:
 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?

2012-03-15, 02:04   #9
LaurV
Romulan Interpreter

"name field"
Jun 2011
Thailand

3·23·149 Posts

Quote:
 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.

2012-03-15, 06:54   #10
Raman
Noodles

"Mr. Tuch"
Dec 2007
Chennai, India

4E916 Posts

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

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

Last fiddled with by Raman on 2012-03-15 at 07:30

2012-03-15, 08:49   #11
Gammatester

Mar 2009

2×19 Posts

Quote:
 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.

Last fiddled with by Gammatester on 2012-03-15 at 08:50

 Similar Threads Thread Thread Starter Forum Replies Last Post a1call Miscellaneous Math 42 2017-02-03 01:29 jasong Miscellaneous Math 5 2016-04-24 03:40 MattcAnderson Puzzles 7 2014-08-17 07:20 kpatsak Math 4 2007-10-29 17:53 m_f_h Puzzles 45 2007-06-15 17:46

All times are UTC. The time now is 22:47.

Mon Feb 6 22:47:25 UTC 2023 up 172 days, 20:15, 1 user, load averages: 1.18, 1.07, 1.06