プロジェクトオイラー
http://projecteuler.net/index.php?section=problems&id=73
前にこの問題を解いた頃は、解ければいいや、ということでずいぶん雑にコードを書いていた。この日は9問解いたようだ。
とにかく何の工夫もなくしらみつぶしに調べていた。
from itertools import imap
from fractions import gcddef count_n(d):
return sum(imap(lambda n: gcd(n, d) == 1,
xrange(d / 3 + 1, (d + 1) / 2)))def solve():
return sum(imap(count_n, xrange(2, N + 1)))N = 12000
print solve()
分母を固定し、分子を変えて、互いに素だったらカウントする。これで40秒近くかかる。
しかし、次のように考えると速くなる。
例えば、分母を60とすると、まず次のような分子が考えられる。
21 22 23 24 25 26 27 28 29
(60 - 1) / 2 - 60 / 3 = 9個ある。
しかし、分子が2の倍数だと約分できるので、これらは取り除く。
22 24 26 28
(60 - 1) / (2 * 2) - 60 / (3 * 2) = 4個ある。
同様に3の倍数、5の倍数も取り除く。
しかし、6の倍数は重複しているので、この分の個数はまた追加する。
このような手順を組み込むと速くなる。
from itertools import imap, ifilterdef make_primes(max_n):
def is_prime(n):
for p in primes:
if p * p > n:
return True
elif n % p == 0:
return False
return True
primes = [ 2 ]
for n in ifilter(is_prime, xrange(3, max_n + 1, 2)):
primes.append(n)
return primesdef factorize(max_n):
b = range(max_n + 1)
a = [ () for n in b ]
for p in primes:
for n in xrange(p, max_n + 1, p):
b[n] /= p
a[n] += (p, )
while b[n] % p == 0:
b[n] /= p
for n in xrange(max_n + 1):
p = b[n]
if p != 1:
a[n] += (p, )
return adef gen_divisors(a, k0 = 0):
yield 1
for k in xrange(k0, len(a)):
d = a[k]
for d2 in gen_divisors(a, k + 1):
yield -d * d2def counter_div(n, d):
if d > 0:
return (n - 1) / (2 * d) - n / (3 * d)
else:
d = -d
return -( (n - 1) / (2 * d) - n / (3 * d))def count_n(d):
counter = 0
return sum(imap(lambda d1: counter_div(d, d1), gen_divisors(f[d])))def solve():
return sum(imap(count_n, xrange(2, N + 1)))N = 12000
primes = make_primes(int(N ** 0.5 + 0.5))
f = factorize(N)
print solve()
こうしたら、0.4秒ほどになった。
このように何も考えないでも解けるけれど、工夫すれば速くなる問題はたくさんあるのではないかと思う。
Forumの最初のほうを見たら、何も考えていないバージョンばかりだった。
