若p是素数且能整除$R(n)$，则显然它也能整除$9R(n)=10^n-1$，即
$$10^{n}\equiv 1\pmod p.$$
显然$p=2$和$p=5$不满足条件，对其他素数p，由费马小定理知$10^{p-1}$模$p$余$1$，因此只需令$a=10^9 \bmod (p-1)$，再检查$10^a \bmod p$是否为1即可知p是否满足条件。重复这一过程直到找到40个素数为止。
由于$a$仍然可能很大，此处采用类似快速幂的分治算法。这40个素数分别是
$$11, 17, 41, 73, 101, 137, 251, 257, 271, 353, 401, 449, 641, 751, 1201, 1409, 1601, 3541, 4001, 4801, 5051, 9091, 10753, 15361, 16001, 19841, 21001, 21401, 24001, 25601, 27961, 37501, 40961, 43201, 60101, 62501, 69857, 76001, 76801, 160001.$$
她们的和是$843296$。
target=1000000
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(2,target) if d[k]==0]
primeSet=set(primeList)
r=set([k for i in range(10) for j in range(10) for k in [(2**i)*(5**j)+1] if k in primeSet and k>5])
def bigMod(k,p):
if k < 10:
return (10**k)%p
else:
return (bigMod(k/2,p)*bigMod(k-k/2,p))%p
i=3
while len(r) < 40:
p=primeList[i]
if p not in r:
k=(10**9)%(p-1)
if bigMod(k,p)==1:
r.add(p)
i+=1
print sum(r)
|
