Project Euler 203

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

Q203.
パスカルの三角形の最初の51列で自乗の因子を含まない異なる自然数の和。

どうやっても解けるような気がする。
素因数分解はあまりやりたくないので、あらかじめ1〜50の素因数分解を行い、その形で組合せの計算を行った。

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 range(3, limit + 1, 2):

        if is_prime(n):

            primes.append(n)
def mul(a, b):

    return [ x + y for x, y in zip(a, b) ]
def div(a, b):

    return [ x - y for x, y in zip(a, b) ]
def value(a):

    n = 1

    for p, e in zip(primes, a):

        if e:

            n *= p ** e

    return n
def is_square_free(a):

    for e in a:

        if e >= 2:

            return False

    return True
def factorial(n):

    return reduce(mul, f[1:n+1])
def C(n, m):

    return div(div(factorial(n), factorial(m)), factorial(n - m))
def factorize(n, i = 0):

    if n == 1:

        return [ 0 ] * len(primes)

    

    for k in range(i, len(primes)):

        p = primes[k]

        if n % p == 0:

            a = [ 0 ] * len(primes)

            a[k] = 1

            return mul(a, factorize(n / p, k))




N = 51

primes = [ 2 ]

make_primes(N)

f = [ 0 ] + [ factorize(n) for n in range(1, N) ]

s = set()

s.add(1)

for n in range(1, N):

    for m in range(1, n / 2 + 1):

        a = C(n, m)

        if is_square_free(a):

            s.add(value(a))

print sum(s)