Project Euler 231

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

Q231.
例えば、120を素因数分解して、120 = 2 × 2 × 2 × 3 × 5 -> 2 + 2 + 2 + 3 + 5 = 14とする。_20000000C_15000000に対して同様の和を求めよ。

n!についてその和を求める。エラトステネスのふるい的に、配列を用意して√20000000までの素数で割っていく。割ったらその分別に用意した変数に足していく。最後に1でない配列の要素は素数だから、そのまま足す。
本当は、20000000!/15000000!についてと5000000!についての和を求めて引き算すればよいが、コードの簡潔さで、階乗についての和の関数を用意するだけにした。

from math import sqrt

from itertools import takewhile
def is_prime(n):

    for p in primes:

        if p * p > n:

            return True

        elif n % p == 0:

            return False

    return True
def make_primes(n):

    m = int(sqrt(n + 0.5))

    for p in xrange(3, m + 1, 2):

        if is_prime(p):

            primes.append(p)

    

    a = [ 0 ] * ((n + 15) / 16)

    for p in primes:

        if p * p > n:

            break

        if p == 2:

            continue

        for k in xrange(3 * p, n, 2 * p):

            a[k>>4] |= 1 << (k & 15)

    

    for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):

        if (a[k>>4] & (1 << (k & 15))) == 0:

            primes.append(k)
def fac_sum(n):

    s = 0

    a = range(n + 1)

    r = int(sqrt(n + 0.5))

    for p in takewhile(lambda p: p <= r, primes):

        p_pow = p

        while p_pow <= n:

            for k in xrange(p_pow, n + 1, p_pow):

                a[k] /= p

                s += p

            p_pow *= p

    

    for e in a:

        if e != 1:

            s += e

    return s
N = 2 * 10 ** 7

M = 15 * 10 ** 6

primes = [ 2 ]

make_primes(int(N ** 0.5 + 0.5))

print fac_sum(N) - fac_sum(M) - fac_sum(N - M)