A slightly revised effort is now up for review as trac #6509: http://trac.sagemath.org/sage_trac/ticket/6509
-Marshall On Jul 10, 2:08 pm, Marshall Hampton <hampto...@gmail.com> wrote: > 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 <sarkar.santanu....@gmail.com> > 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 sage-support@googlegroups.com To unsubscribe from this group, send email to sage-support-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-support URLs: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
