Project Euler 154

プロジェクトオイラー
http://projecteuler.net/index.php

Q154.
(x + y + z)200000の係数で、それが1012の倍数である個数

n = m + k + lとして、n!/m!k!l!が係数となる。それが2と5を12個以上含めばよい。
5の個数は、[n/5] + [n/52] + ... - ([k/5] + [k/52] + ... + [l/5] + [l/52] + ... + [m/5] + [m/52] + ...)となる。そのあとあれこれ考えたが、結局素直にやった。しかし遅い。が、k,l,mは正三角形を成していて、その対称性より、だいたい1/6計算すればよい。



import time

def P(n, p):
s = 0
while n:
n /= p
s += n
return s

def P_min(n, p):
if p == 2:
s = P2[n]
elif p == 5:
s = P5[n]
else:
s = P(n, p)

while n >= p:
if n % p != p - 1:
break
n /= p
while n >= p:
s -= 1
n /= p
return s

t = time.clock()
N = 200000
M = 12

P2 = [ P(n, 2) for n in xrange(N + 1) ]
P5 = [ P(n, 5) for n in xrange(N + 1) ]

counter = 0
n2 = P2[N]
n5 = P5[N]
for k in xrange((N + 2) / 3, N + 1):
if (k & 255) == 255:
print k
k2 = P2[k]
k5 = P5[k]
lm = N - k
lm_min2 = P_min(lm, 2)
lm_min5 = P_min(lm, 5)
if n5 - k5 - lm_min5 < M:
continue
elif n5 - k5 - P5[lm] >= M:
counter += lm + 1
continue

for l in xrange((lm + 1) / 2, min(k, lm) + 1):
m = lm - l
lm5 = P5[l] + P5[m]
if n5 - k5 - lm5 < M:
continue
lm2 = P2[l] + P2[m]
if n2 - k2 - lm2 < M:
continue
if k == l or l == m:
counter += 3
else:
counter += 6

print counter
print time.clock() - t