满足要求的素数共有16个，分别是5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149。
对每个这样的素数p，找出满足$a_p^2+b_p^2=p$的解$a_p,b_p$，本题中这样的解是唯一的。
利用公式
$$(a^2+b^2)(c^2+d^2)=(ac-bd)^2+(ad+bc)^2=(ac+bd)^2+(ad-bc)^2,$$
即可使用深度优先搜索枚举这些素数的乘积并生成不同的解。最终结果是$2032447591196869022$。
注：计算量较大，代码中打印了进度信息。
primes=[5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149]
squares=[i*i for i in range(1,13)]
ap,bp=[],[]
for p in primes:
for j,x in enumerate(squares):
y=p-x
if y in squares and 0 < x <= y:
k=squares.index(y)
ap.append(j+1)
bp.append(k+1)
q=[(p,{(ap[i],bp[i])}) for i,p in enumerate(primes)]
checked=set(primes)
r=sum(ap)
k=16
while q:
u,s=q.pop()
for i in range(16):
if u%primes[i]>0:
c,d,v=ap[i],bp[i],u*primes[i]
if v not in checked:
checked.add(v)
solutions=set()
for a,b in s:
for x,y in ((abs(a*c-b*d),abs(a*d+b*c)),(abs(a*c+b*d),abs(a*d-b*c))):
if x>y:
x,y=y,x
if (x,y) not in solutions:
solutions.add((x,y))
r+=x
q.append((v,solutions))
k+=1
if k%655==0:
print k/655,"percent completed"
print r
|
