依次考虑能分别整除$p, p^2, \ldots$的数个数就可以看出$n!$中$p$的幂是
$$\lfloor\frac{n}{p}\rfloor+\ldots+\lfloor\frac{n}{p^k}\rfloor$$
其中k满足$p^k\le n < p^{k+1}$。先用筛法生成小于$10^6$的所有素数,以及每个数的素因子列表。
对每个i,用上述数据找出$i!$中幂次最高的素数及其幂,将幂次乘以1234567890后用二分查找确定满足要求的$N(i)$,最终结果是$278157919195482643$。
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]
print "prime list generated"
factors=[[] for i in range(target+1)]
i=0
while i < len(primeList) and primeList[i] <= target:
p=primeList[i]
while p < n:
factors[p].append(primeList[i])
p+=primeList[i]
i+=1
print "factor list generated"
def getPowers(n,p):
t=0
while n%p==0:
n/=p
t+=1
return t
def getPowersFact(n,p):
t,m=0,n
while m>0:
m/=p
t+=m
return t
def getN(n,p):
threshold=1234567890*getPowersFact(n,p)
m=(p-1)*threshold
m+=p-m%p
t=0
bound=threshold-getPowersFact(m,p)
while t < bound:
m+=p
t+=getPowers(m,p)
return m
i=11
q=5
r=9876543150
for i in range(11,target+1):
x=getN(i,q)
a,b=max((getN(i,p),p) for p in factors[i])
if a>x:
x,q=a,b
r=(r+x)%1000000000000000000
print r
|