On Fri, Jul 10, 2009 at 10:01 AM, Marshall Hampton<[email protected]> wrote: > > 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))]
There is a function two_squares in Sage already: sage: two_squares(920834098234) Yours above miht be better than it though... > > 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 > > > -- William Stein Associate Professor of Mathematics University of Washington http://wstein.org --~--~---------~--~----~------------~-------~--~----~ 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 -~----------~----~----~----~------~----~------~--~---
