Sum of Two Squares

[Eddie W. Shore]

Given a positive integer n, can we find two non-negative integers x and y such that:

n = x^2 + y^2

The program presented here is the use of iterations to find all possible pairs which fit n = x^2 + y^2. Some integers do not have representations, others have more than one. The program will show all possible combinations.

HP Prime Program SUM2SQ

Code:


EXPORT SUM2SQ(n)
BEGIN
// EWS 2019-07-21
// breaking n into a sum of 2 squares
LOCAL r,j,k,l;
// we can more than 1 representation
r:=IP((n/2)^0.5);
l:={};
FOR j FROM 0 TO r DO
k:=(n-j^2)^0.5;
IF FP(k)==0 THEN
l:=CONCAT(l,
{STRING(j)+"^2 + "+
STRING(k)+"^2 = "+
STRING(n)});
END;
END;

RETURN l;
END;

Blog link: http://edspi31415.blogspot.com/2019/08/h...f-two.html

[klesl]

Very interesting :-)
Another interesting topic in number theroy is Bézout's_identity
https://en.wikipedia.org/wiki/Bézout's_identity
Can you write a program for this also?

[Eddie W. Shore]

I'll check it out and maybe give it a shot.

[Eddie W. Shore]

(08-17-2019 07:05 PM)klesl Wrote:  Very interesting :-)
Another interesting topic in number theroy is Bézout's_identity
https://en.wikipedia.org/wiki/Bézout's_identity
Can you write a program for this also?

Here is the Bezout program: this program returns the gcd of two integers a and b, and all of the possible integer solutions that fit a*x + b*y = d, with in a given range [-m, m]:

Code:

EXPORT BEZOUT(a,b,m)
BEGIN
// integers a, b, range -m to m
// 2019-08-19 EWS
LOCAL d,l,k,x,y;
d:=gcd(a,b);
l:={};
FOR k FROM −m TO m DO
x:=k;
y:=(d-a*x)/b;
IF FP(y)==0 THEN
l:=CONCAT(l,(x,y));
END;
END;
RETURN {d,l};
END;

Output: { gcd, list of (x,y) pairs }

Example: a = 18, b = 45, m = 20, find all solutions within the range [-20, 20]:
BEZOUT(18, 45, 20):

{9, {(-17, 7), (-12,5), (-7,3), (-2,1), (3,-1), (8,-3), (13,-5), (18,-7)}}

[Albert Chan]

(08-17-2019 01:35 PM)Eddie W. Shore Wrote:  Given a positive integer n, can we find two non-negative integers x and y such that:

n = x^2 + y^2

If we have 2 pairs,, we can factor n

HOME> SUM2SQ(12349)

{"30^2 + 107^2 = 12349" , "75^2 + 82^2 = 12349"}

gcd(12349, (75-30)^2 + (107-82)^2)
= gcd(12349, 45^2 + 25^2)         // gcd(12349, 5) = 1, pull 5^2 out
= gcd(12349, 9^2 + 5^2 = 106)   // gcd(12349, 2) = 1, pull 2 out
= gcd(12349, 106/2 = 53)            // 12349 = 53 * 233
= 53

Proof:
n = a^2 + b^2 = c^2 + d^2

(a^2 - c^2) = (d^2 - b^2)
(a-c)*(a+c) = (d-b)*(d+b)            // Let f=a-c, f'=a+c, g=d-b, g'=d+b
f f' = g g'

(f^2 + g^2) * (f'^2 + g^2)
= (f f')^2 + g^2*(f^2 + f'^2) + g^4
= (g g')^2 + g^2*(f^2 + f'^2 + g^2)
= g^2 * (f^2 + f'^2 + g^2 + g'^2)
= g^2 * 4n

-> n = (f^2 + g^2) * (f'^2 + g^2) / (4 g^2)

To keep (f^2 + g^2) small, it is better if (a,b) and (c,d) are in same sort order.
Ref: book "Dead Reckoning: calculating without instruments", p64

[Albert Chan]

(09-09-2021 11:12 PM)Albert Chan Wrote:  n = a^2 + b^2 = c^2 + d^2

(a^2 - c^2) = (d^2 - b^2)
(a-c)*(a+c) = (d-b)*(d+b)            // Let f=a-c, f'=a+c, g=d-b, g'=d+b
f f' = g g'
...
-> n = (f^2 + g^2) * (f'^2 + g^2) / (4 g^2)

Extend this to n = a^2 + k*b^2 = c^2 + k*d^2
Do the math (k=1 should degenerate to above quoted identities)

f f' = k g g'

n = (f^2 + k*g^2) * (f'^2 + k*g^2) / (4k g^2)
   = (2n - 2*(a*c+k*b*d)) * (2n - 2*(a*c-k*b*d)) / (4k g^2)
   = ((a*c + k*b*d) - n) * ((a*c - k*b*d) - n) / (k*g^2)

Numerator must be multiples of n:

(a*c + k*b*d) * (a*c - k*b*d) ≡ 0 (mod n)           // multiply d^2, k*d^2 ≡ -c^2 (mod n)
(a*c*d-b*c^2) * (a*c*d+b*c^2) ≡ 0 (mod n)       // divide by c^2 ≠ 0 (mod n)
(a*d - b*c) * (a*d + b*c) ≡ 0 (mod n/gcd(n,c^2))

gcd(n, a*d ± b*c) produce factor of n

Redo previous example (k=1):

12349 = 30^2 + 107^2 = 75^2 + 82^2

gcd(12349, 30*82-107*75 = -5565) = gcd(5565, 12349-5565 = 6784 = 2^7*53) = 53
gcd(12349, 30*82+107*75 = 10485) = gcd(10485, 12349-10485 = 1864 = 2^3*233) = 233

Now, make k = 15^2 = 225:

12349 = 107^2 + 225*2^2 = 82^2 + 225*5^2

gcd(12349, 107*5-82*2 = 371 = 7*53) = 53
gcd(12349, 107*5+82*2 = 699 = 3*233) = 233

Another example, from "Dead Reckoning", p65, (k = 22)

1457 = 37^2 + 22*2^2 = 7^2 + 22*8^2

gcd(1457, 37*8-2*7 = 282 = 2*3*47) = 47
gcd(1457, 37*8+2*7 = 310 = 2*5*31) = 31

[Albert Chan]

(09-10-2021 02:55 PM)Albert Chan Wrote:  Another example, from "Dead Reckoning", p65, (k = 22)

1457 = 37^2 + 22*2^2 = 7^2 + 22*8^2

gcd(1457, 37*8-2*7 = 282 = 2*3*47) = 47
gcd(1457, 37*8+2*7 = 310 = 2*5*31) = 31

Another way, by forcing k=1 ... radical will disappear.

1457 = 37^2 + (2√22)^2 = 7^2 + (8√22)^2

gcd(1457, 30^2+(6√22)^2 = 1692 = 1457+5*47) = 47

This work even for sum of square of radicals.

1457 = 3*2^2 + 5*17^2 = (2√3)^2 + (17√5)^2
1457 = 3*22^2 + 5*1^2 = (22√3)^2 + (√5)^2

gcd(1457, (20√3)^2+(16√5)^2 = 2480 = 2^4*5*31) = 31

Interestingly, gcd(n, a*d ± b*c) seems to work as well.

gcd(1457, 2*1-17*22 = -372 = -2^2*3*31) = 31
gcd(1457, 2*1+17*22 = 376 = 2^3*47) = 47