Project Euler 88

PE 数学 プログラミング

Problem 88
自然数Nが少なくとも2つの自然数の集合{a1, a2, ... , ak}の和と積で書けるとき、これを積和数と呼ぶ：N = a1 + a2 + ... + ak = a1 × a2 × ... × ak。
例えば、6 = 1 + 2 + 3 = 1 × 2 × 3。
与えられた大きさkとの集合に対して、この性質を持つ最小のNを最小積和数と呼ぶ。大きさの集合、k = 2, 3, 4, 5, 6の最小積和数は次のようになる。
k=2: 4 = 2 × 2 = 2 + 2
k=3: 6 = 1 × 2 × 3 = 1 + 2 + 3
k=4: 8 = 1 × 1 × 2 × 4 = 1 + 1 + 2 + 4
k=5: 8 = 1 × 1 × 2 × 2 × 2 = 1 + 1 + 2 + 2 + 2
k=6: 12 = 1 × 1 × 1 × 1 × 2 × 6 = 1 + 1 + 1 + 1 + 2 + 6
ゆえに、2≤k≤6に対して、全ての最小積和数の和は4+6+8+12 = 30。ここで、8は和の中で一度しかカウントしていないことに注意。
ちなみに、2≤k≤12に対する最小積和数の完全な集合は{4, 6, 8, 12, 15, 16}なので、和は61となる。
2≤k≤に対する最小積和数の総和をいくつか。
http://projecteuler.net/index.php?section=problems&id=88

まず、素数は積和数にはなりえません。
合成数の場合、例えば6を考えると、

2 × 3 - (2 + 3) = 1

で、和が足りない分は1で補い、

1 × 2 × 3 = 1 + 2 + 3

k = 3となります。
順番に素因数分解して、積を出します。このときのkを出して、

n 積 k
4 2 × 2 2 ○
6 2 × 3 3 ○
8 2 × 4 4 ○
8 2 × 2 × 2 5 ○
9 3 × 3 5
10 2 × 5 5
12 2 × 6 6 ○
12 2 × 2 × 3 8 ○
12 3 × 4 7 ○

○がついたところがkに対する最小の積和数です。すべてのkについてnを求めます（Code1）。

これよりある数以下の積を出すほうが速いです。この方法だと素数も素因数分解もする必要がありません。ただし、いくつまでの積を出していいのかわからないので、倍々にしていきます（Code2）。

# Code1
from itertools import count, imap, ifilter, takewhile

def head(a):
    for e in a:
        return e

def is_prime(n):
    return all(n % p != 0 for p in takewhile(lambda p: p * p <= n, primes))

def gen_primes():
    for k in count():
        if k == len(primes):
            primes.append(head(ifilter(is_prime, count(primes[-1] + 1))))
        yield primes[k]

def div_pow(n, p):
    e = 0
    while n % p == 0:
        n /= p
        e += 1
    return n, e

def factorize(n):
    if n == 1:
        return [ ]
    
    for p in takewhile(lambda p: p * p <= n, gen_primes()):
        n, e = div_pow(n, p)
        if e > 0:
            return [(p, e)] + factorize(n)
    
    return [(n, 1)]

def value(f):
    return reduce(lambda x, y: x * y[0] ** y[1], f, 1)

def gen_divisors(f, k = 0):
    def add(f, p, e):
        if e == 0:
            return f
        else:
            return ((p, e),) + f
    
    if k == len(f):
        yield (), ()
    else:
        p, e0 = f[k]
        for f1, f2 in gen_divisors(f, k + 1):
            for e in range(e0 + 1):
                yield add(f1, p, e), add(f2, p, e0 - e)

def gen_product(f, limit = 2):
    for f1, f2 in gen_divisors(f):
        d1 = value(f1)
        d2 = value(f2)
        if limit <= d1 <= d2:
            yield d1 + d2, 2
            for s, l in gen_product(f2, d1):
                yield d1 + s, l + 1

def calc_k(N):
    def is_full():
        return all(a[k] != 0 for k in xrange(N, 1, -1))
    
    a = [ 0 ] * (N + 1)
    for n in takewhile(lambda n: not is_full(), count(4)):
        f = factorize(n)
        for k in (n + l - s for s, l in gen_product(f)):
            if k <= N and a[k] == 0:
                a[k] = n
    return a

N = 12 * 10 ** 3
primes = [ 2 ]
print sum(set(calc_k(N)))

# Code2
from itertools import count, imap, ifilter, takewhile

def head(a):
    for e in a:
        return e

def gen_sum_product(begin, end, s = 0, p = 1, limit = 2):
    for n in xrange(limit, end):
        if s:
            yield s + n - 1, p * n
        for s1, p1 in gen_sum_product((begin - 1) / n + 1,
                        (end - 1) / n + 1, s + n - 1, p * n, n):
            yield s1, p1

def calc_k(N):
    def is_full():
        return all(a[k] != 0 for k in xrange(N, 1, -1))
    
    a = [ 0 ] * (N + 1)
    begin, end = 1, 2
    while not is_full():
        for s, n in gen_sum_product(begin, end):
            k = n - s
            if k <= N and (a[k] == 0 or a[k] > n):
                a[k] = n
        begin, end = end, end * 2
    
    return a

N = 12 * 10 ** 3
primes = [ 2 ]
print sum(set(calc_k(N)))