OEIS A216451中给出了本问题的同余性条件，不过由于本题数据范围非常大，很难有效筛选范围内的素数，所以似乎并不存在能有效利用此条件的手段（见注2）。
以下代码使用的是最朴素的思想：枚举a,b以生成所有形如$a^2+db^2$的数并在$d=1,2,3,7$中取交集。
由于内存有限，所以需要将20亿以内的数按每5百万一组分成400组依次筛查，并用列表f记忆对每个a，b所枚举到的位置。
最终结果是11325263。
注：由于计算量很大（大概需要运行一小时），代码中打印了进度信息。
注2：找出形如$x^2+dy^2$素数是一个经典问题，可以参考该页面以及其中提到的参考书《Primes of the form $x^2+ny^2$》。
对于具体的d，可以从该索引页面找到更进一步的信息。
假设n是一个满足条件的数，将其写成$n=n_0m^2$，其中$n_0$是无平方因子数。
设p是$n_0$的一个素因子，则p能整除$x^2+dy^2$形式的数，此处$d=1,2,3,7$，此时有经典结论：-d是p的二次剩余。
对$d=1,2,3,7$的情况，上述情况还可以进一步简化为如下的同余性条件：
$d=1$： $p=2$或p模4余1。
$d=2$： $p=2$或p模8余1,3。
$d=3$： $p=3$或p模3余1。
$d=7$： p模28余1, 7, 9, 11, 15, 23, 25。
综合上述条件可知p需要同时满足模8余1（此时模4也余1，另外模8余3与模4余1冲突）、模3余1、模28余1, 9, 25（余7,11,15,23的情况与模4余1的情况冲突）。
取$c=\mathrm{lcm}(3,8,28)=168$，检验可知p需要是模168余1,25,121的素数。
因此需要先找出这样的素数，并组合成满足条件的无平方因子数，再与满足条件的完全平方数组合，简单试算可见复杂度并不显著小于上文中的朴素方法。
target=2000000000
block=5000000
bound=44722
squares=[i*i for i in range(bound)]
f1=[1]*bound
f2=[1]*bound
f3=[1]*bound
f7=[1]*bound
r=0
for left in range(0,target,block):
flags=[0]*block
right=left+block
i=1
while i < bound and squares[i] < right:
for d,f,k in [(1,f1,1),(2,f2,2),(3,f3,4),(7,f7,8)]:
j=f[i]
while j < bound and squares[i]+d*squares[j] < right:
flags[squares[i]+d*squares[j]-left]|=k
j+=1
f[i]=j
i+=1
r+=flags.count(15)
print right/block, "out of 400 done,", r, "collected"
print r
|
