I tried to do this for fun, but my code doesn't quite work. There is
a java applet with quite readable source code available at:
http://www.alpertron.com.ar/FSQUARES.HTM
I was trying to implement the ideas at: http://www.schorn.ch/howto.html
but I am out of time for now, maybe someone else can fix this:
def sum_of_two_squares(p):
if mod(p,8)==5:
b = 2
else:
pnum = 2
b = nth_prime(pnum)
q = mod(b^((p-1)/2),p)
while q == 1:
pnum += 1
b = nth_prime(pnum)
q = mod(b^((p-1)/2),p)
b = mod(b^((p-1)/4),p)
a = p
while b^2 > p:
a,b = b, mod(a,b)
return [Integer(b), Integer(mod(a,b))]
def sq4(n):
if squarefree_part(Integer(n)) == 1:
return [0,0,0,sqrt(n)]
m = n
v = 0
while mod(m,4)==0:
v = v +1
m = m/4
if mod(m,8) == 7:
d = 1
m = m - 1
else:
d = 0
if mod(m,8)==3:
x = int(sqrt(m))
if mod(x,2) == 0:
x = x - 1
p = (m-x^2)/2
while not is_prime(p):
x = x - 2
p = (m-x^2)/2
y,z = sum_of_two_squares(p)
return [2^v*q for q in [d,x,y+z,abs(y-z)]]
x = int(sqrt(m))
p = m - x^2
if p == 1:
return[2^v*q for q in [d,0,x,1]]
while not is_prime(p):
x = x - 1
p = m - x^2
y,z = sum_of_two_squares(p)
return [2^v*q for q in [d,x,y,z]]
which works for most numbers, but not all.
-M. Hampton
On Jul 10, 8:30 am, Santanu Sarkar <[email protected]>
wrote:
> Is the algorithm for sum of four square i,e every positive integer
> can be expressed as the sum of four square is implemented in
> Sage
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
