Sums of two squares by rolf

f=lambda u,n:f(f4(u),n-1)if n else u
def E(a,b,n):
 if a<n**.5:return(b,a)
 if b:return E(b,a%b,n)
 return a
def M(x):
 h=x[0];o=x[1:] 
 while o:f=o[0];t=o[1:];h=set(sum(([f1(a,b),f2(a,b)]for a in h for b in f),[]));o=t
 return h  
def f0(a,b):return tuple(sorted((abs(a),abs(b))))
def f1(u,v):a,b=u;c,d=v;a,b=a*c+b*d,a*d-b*c;return f0(a,b)
def f2(u,v):a,b=u;c,d=v;a,b=a*c-b*d,a*d+b*c;return f0(a,b)
def f3(u):a,b=u;a,b=a**2-b**2,2*a*b;return f0(a,b)
def f4(u):a,b=u;a,b=a+b,a-b;return f0(a,b)
def f5(u):a,b=u;a,b=0,a**2+b**2;return f0(a,b)
def U(p):
 p,e=p;q=2
 while any(x**2%q==p%q for x in range(1,q+1)):q+=1
 q=E(p,pow(q,(p-1)/4,p),p)
 if e==1:return set([q])
 a,b=m=f3(q),f5(q)
 if e==2:return set(m)
 if e==3:return set(f1(w,q)for w in m)|set(f2(w,q)for w in m)
 m=set([f1(a,b),f2(a,b),f3(a),f3(b)])
 while e>4:m=set(f1(w,q)for w in m)|set(f2(w,q)for w in m);e-=1
 return m
def V(p):p,e=p;return set([(0,p**(e/2))])
while 1:
 n=input();p=input();S=0
 T=p[0][1]if p[0][0]==2 else 0
 g=filter(lambda x:(x[0]-1)%4==0,p) 
 h=filter(lambda x:(x[0]-3)%4==0,p)
 if all(x[1]%2==0 for x in h):
  if g or h:
   o=[U(p)for p in g]+[V(p)for p in h]
   S=[f(x,T)for x in M(o)]
  else:   
   t=(2**T/2)**.5
   S=[0,set([(int(t),int(t))])][t==int(t)]
 if S:print" = ".join("%d^2 + %d^2"%(x,y)for x,y in sorted(S))
 else:print

Note that non-ascii characters in the above source code will be escaped (such as \x9f).

To protect the system from spam, please input your favorite sport (hint: I believe its name must start with 'g', case insensitive)

download

return to the top page