——–求sigma{S(N)}，其中N分解因素{zg}幂为1，N只能被4k+1形式的素数整除，S(N)为N = a^2+b^2(0<=a<=b)所有a的和
——–所有符合条件的素数递归枚举得出所有N，N = a^2+b^2，M = c^2+d^2，NM = (ac)^2+(ab)^2+(bc)^2+(bd)^2=(ac + bd)^2+(ad – bc)^2=(ac – bd)^2 + (ad + bc)^2. 所以每乘一个素数的到N = an^2 + bn^2, 所以an = min(ac+bd, abs(ad-bc))或 an = min(abs(ac-bd), ad+bc).

import math
MAXN = 150
prime = [0 for x in range( MAXN )]
np = 0
def primelist():
? ? global prime, np
? ? for i in range( MAXN ):
? ? ? ? prime[ i ] = 1
? ? for i in range(2, MAXN):
? ? ? ? if( prime[ i ] ):
? ? ? ? ? ? for j in range( i*i, MAXN, i ):
? ? ? ? ? ? ? ? ? ? prime[ j ] = 0
? ? j = 0
? ? for i in range(2, MAXN):
? ? ? ? if( prime[ i ] and (i - 1)%4 == 0 ?):
? ? ? ? ? ? prime[ j ] = i
? ? ? ? ? ? j += 1
? ? np = j

primelist()
a = [ 0 for x in range( np )]
b = [ 0 for x in range( np )]

for i in range( np ):
? ? temp = int(math.sqrt(prime[ i ]/2)) + 1
? ? for j in range( temp ):
? ? ? ? rest = prime[ i ] - j*j
? ? ? ? sq = int(math.sqrt( rest ))
? ? ? ? if( sq*sq == rest ):
? ? ? ? ? ? a[ i ] = j
? ? ? ? ? ? b[ i ] = sq

def imin( x, y):
? ? if(x < y ):
? ? ? ? return x
? ? return y
def imax( x, y):
? ? if( x > y ):
? ? ? ? return x
? ? return y

sum = 0
def getsum(p, pro, ap, bp):
? ? global prime, np, sum, a, b
? ? sum += ap
? ? for i in range( p + 1, np ):
? ? ? ? a1temp = abs(ap*a[ i ] + bp*b[ i ])
? ? ? ? b1temp = abs(ap*b[ i ] - bp*a[ i ])
? ? ? ? a2temp = abs(ap*a[ i ] - bp*b[ i ])
? ? ? ? b2temp = abs(ap*b[ i ] + bp*a[ i ])
? ? ? ? a1next = imin(a1temp, b1temp)
? ? ? ? b1next = imax(a1temp, b1temp)
? ? ? ? a2next = imin(a2temp, b2temp)
? ? ? ? b2next = imax(a2temp, b2temp)
? ? ? ? getsum( i, pro * prime[ i ], a1next, b1next )
? ? ? ? if( a1next != a2next ):
? ? ? ? ? ? getsum( i, pro * prime[ i ], a2next, b2next )

for i in range( np ):
? ? getsum( i, prime[ i ], a[ i ], b[ i ] )

print sum