I have tested this up to 10000, seems to work:
def brute_force_s4(n):
'''
Brute force search for decomposition into a sum of four squares,
for cases that the main algorithm fails to handle.
'''
for i in range(0,int(sqrt(n))+1):
for j in range(i,int(sqrt(n-i^2))+1):
for k in range(j, int(sqrt(n-i^2-j^2))+1):
l = int(sqrt(n-i^2-j^2-k^2))
if n-i^2-j^2-k^2-l^2==0:
return [i,j,k,l]
def sq4(n):
'''
Computes the decomposition into the sum of four squares,
using an algorithm described by Peter Schorn at:
http://www.schorn.ch/howto.html.
'''
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
if x < 0:
# fall back to brute force
m = m + d
return [2^v*q for q in brute_force_s4(m)]
y,z = 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
if x < 0:
# fall back to brute force
m = m + d
return [2^v*q for q in brute_force_s4(m)]
y,z = two_squares(p)
return [2^v*q for q in [d,x,y,z]]
Cheers,
Marshall Hampton
On Jul 10, 9: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
-~----------~----~----~----~------~----~------~--~---
