Project Euler 193

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

Q193.
平方数で割り切れない2⁵⁰未満の自然数の個数。

割り切れるほうの個数を考える。n²で割り切れる自然数の集合の個数をS_nとすると、なんらかの平方数で割り切れる自然数集合の個数Sは、和集合の公式を使って、

S = S₂ + S₃ + S₅ + ... - S_2・3 - S_2・5 - ... + S_2・3・5 + ...

2²⁵までのこのような複数の同じ素数を使わない自然数を列挙する。その際、奇数個の素数を使う自然数は負に変えておくと、あとの計算が便利になる。

def sieve(n):

    a = [ False, False ] + [ True ] * (n - 1)

    p = 2

    while p * p <= n:

        for q in xrange(p * 2, n + 1, p):

            a[q] = False

        p = a.index(True, p + 1)

    

    return filter(lambda i: a[i], xrange(n + 1))
def gen_prime_combinations(n, k = 0):

    for l in xrange(k, len(primes)):

        p = primes[l]

        if p > n:

            break

        yield -p

        m = n / p

        if m > 1:

            for a in gen_prime_combinations(m, l + 1):

                yield -p * a
def mul(x, y): return x * y
M = 2 ** 25

N = M * M

primes = sieve(M)

s = N

for a in gen_prime_combinations(M):

    if a > 0:

        s += N / (a * a)

    else:

        s -= N / (a * a)

print s