Project Euler 53

Problem 53
（前略）1 ≤ n ≤ 100 で_nC_rの値が100万を超えるものは区別せずにいくつあるか。
http://projecteuler.net/index.php?section=problems&id=53

_nC_rの計算は、

_nC_r = _nC_r-1 * (n - r + 1) / r

という漸化式で計算します。

from itertools import imap
def C(n, r):

    return 1 if r == 0 else C(n, r - 1) * (n - r + 1) / r
def gen_combination():

    for n in xrange(N + 1):

        for r in xrange(n + 1):

            yield C(n, r)
N = 100

M = 10 ** 6

print sum(imap(lambda n: n > M, gen_combination()))

しかし、これだと例えば₇C₂は₇C₃、₇C₄などでも計算することになりムダです。次のようにすればこのムダはなくなり、一気に速くなります。

from itertools import imap
def gen_combination():

    for n in xrange(N + 1):

        C = 1

        for r in xrange(n + 1):

            yield C

            C = C * (n - r) / (r + 1)
N = 10 ** 2

M = N ** 3

print sum(imap(lambda n: n > M, gen_combination()))

実はこんなにたくさん計算する必要はありません。パスカルの三角形で考えて、上下や左右の大小関係はわかっているので、100万の境界だけ分かれば数えることができます。具体的には、まず₁₀₀C₀、₁₀₀C₁、と計算していき、₁₀₀C_rがはじめて100万を超えるrを見つけます。そうすると、このr以上で中央を挟んで反対側まで100万を超える範囲なので、(n-2r + 1)個が100万を超えます。次に、₉₉C_rを計算し、100万以下ならrを大きくし、そうでなければ小さくして境界を探します。境界が見つからなくなったらそこで終わりです。

from itertools import imap
def search_up(n, r0, C):

    for r in xrange(r0 + 1, n / 2 + 1):

        C = C * (n - r + 1) / r

        if C > M:

            return r, C

    return -1, 0
def search_down(n, r0, C0):

    for r in xrange(r0 - 1, -1, -1):

        C = C0 * (r + 1) / (n - r)

        if C <= M:

            return r + 1, C0

        C0 = C
# minimum r s.t. nCr > M

def min_r(n, C, r0):

    if C == 1:

        return search_up(N, 0, 1)

    else:

        C = C * (n - r0 + 1) / (n + 1)

        if C <= M:

            return search_up(n, r0, C)

        else:

            return search_down(n, r0, C)
N = 10 ** 2

M = N ** 3

counter = 0

C = 1

r = 0

for n in xrange(N, -1, -1):

    r, C = min_r(n, C, r)

    if r == -1:

        break

    counter += n + 1 - 2 * r

print counter