Project Euler 174

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

Q174.
前の問題で、同じタイルの枚数を並べ替えて別の形状にできることがある。これが10通り以下のタイルの枚数は100万までにいくつあるか。

タイルの枚数は、4m(n-m) (m<n-m)だから、m(n-m)の約数が21個以下なら10通り以下。約数の個数は素因数分解の形で決まる。21以下になる形は、p₁p₂p₃p₄など、数が限られるので、それを分割数と同じように生成し、その個々の形について、100万以下になる個数を求める。その際に、素数の個数関数π(x)の逆関数のようなものを作っておくと、計算が速くなる。その関数には、π(x) 〜 x / log xを使う。

from math import log
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(limit):

    for n in xrange(3, limit, 2):

        if is_prime(n):

            primes.append(n)
def pi_inverse(n):

    if n == 1:

        return 0

    i_max = len(primes) - 1

    ln_n = log(float(n))

    i1 = min(i_max, int(n / ln_n))

    while True:

        n1 = primes[i1]

        if n1 == n:

            return i1

        elif abs(n1 - n) > 100:

            i1 -= int((n1 - n) / ln_n)

        elif n1 > n:

            while n1 > n:

                i1 -= 1

                n1 = primes[i1]

            return i1

        else:

            while n1 <= n:

                i1 += 1

                n1 = primes[i1]

            return i1 - 1
def gen_prime_form(n, l = -1):

    if l == -1:

        l = n - 1

    elif l >= n:

        l = n - 1

    

    if n > 1:

        for k in range(1, n):

            yield [ k ]

            for a in gen_prime_form(n / (k + 1), k):

                yield [ k ] + a
def f(a, n, i0 = 0):

    m = int((n + 0.5) ** (1. / sum(a)))

    if len(a) == 1:

        if m >= primes[i0]:

            return pi_inverse(m) - i0 + 1

        else:

            return 0

    else:

        s = 0

        for i1 in xrange(i0, pi_inverse(m) + 1):

            s += f(a[1:], n / primes[i1] ** a[0], i1 + 1)

        return s
N = 10 ** 6 / 4

M = 10
primes = [ 2 ]

make_primes(N + 1000)

print sum(map(lambda a: f(a, N), gen_prime_form(M * 2 + 1)))