本题的思路很直接，先生成素数列表。对每个素数p，用Lucas定理计算出$\binom{n}{k}$模p的余数，再用中国剩余定理将其组合起来即可，最终结果是$162619462356610313$。
import itertools
target=5000
d=[0]*target
n=2
while n < target:
for i in range(n*n,target,n):
d[i]=d[i]+1
n=n+1
while n < target and d[n]>0:
n=n+1
primeList=[k for k in range(1000,target) if d[k]==0]
a={}
f=[0]*(target+1)
inverses={}
C=lambda n,k,p:0 if k>n else ((f[n]*inverses[p][f[k]])%p)*(inverses[p][f[n-k]]%p)
lucas=lambda n,k,p:1 if n==0 else (C(n%p,k%p,p)*lucas(n/p,k/p,p))%p
CRT=lambda a,b,c,p,q,r:(a*inverses[p][(q*r)%p]*q*r+b*inverses[q][(p*r)%q]*p*r+c*inverses[r][(p*q)%r]*p*q)%(p*q*r)
for p in primeList:
inverses[p]=[0]*p
inverses[p][1]=1
for i in range(2,p):
inverses[p][i]=((p-p/i)*inverses[p][p%i])%p
f[0]=1
for i in range(1,p):
f[i]=(f[i-1]*i)%p
a[p]=lucas(10**18,10**9,p)
print sum(CRT(a[p],a[q],a[r],p,q,r) for p,q,r in itertools.combinations(primeList,3) if a[p] or a[q] or a[r])
|
