Project Euler 77

Problem 77
10は素数の和に5つの方法で書ける。
7 + 3
5 + 5
5 + 3 + 2
3 + 3 + 2 + 2
2 + 2 + 2 + 2 + 2
素数の和で書く方法が5000を超える最初の値はいくつか。
http://projecteuler.net/index.php?section=problems&id=77

分割数の素数版とも言うべきものです。
前回のCode1の方法を真似れば比較的簡単です（Code1）。

前回のCode2の方法を考えてみましょう。この方法は母関数の係数を何項まで求めるかがわかっていないとつらいです。それを倍、倍としていく方法もあることはありますが（Code2）。
Haskellなら遅延評価で何項まで求めるかわかっていなくても大丈夫です。遅延評価というのは実際に使うときに初めて値を求めるということです。この仕組みをシミュレートしてみましょう。Haskellで書いて、

a = start:[ f e | e <- a ]

これをシミュレートします。
fは要するに次の項を出す関数です。例えばfが1を足した数を返す関数とすると、

a = 1:[ n + 1 | n <- a ]    -- [1, 2, 3, … ]

初期のリストをa、次を出す関数をnext、初期リスト以外に次を出すのに必要な引数をまとめたものをargsとします。上のリストならa = [ 1 ]、argsは適当にします。リストを成長させていきます。

from itertools import count, takewhile

def gen(x):
    a, next, args = x
    for k in count():
        if k < len(a):
            yield a[k]
        else:
            e = next(a, args)
            a.append(e)
            yield e

a = [ 1 ]
for n in takewhile(lambda n: n < 3, gen((a, lambda b, d: b[-1] + 1, 0))):
    print n,
for n in takewhile(lambda n: n < 5, gen((a, lambda b, d: b[-1] + 1, 0))):
    print n,

これで、

1 2 1 2 3 4

と出力されます。
このgenを利用して書いて動かしてみると、再帰が深すぎると怒られてしまいました。
再帰の深さはデフォルトで1000までのようです。これを次のように10000にすると、

sys.setrecursionlimit(10000)

結果が出力されました（Code3）。最初のコードより1桁遅いです。

# Code1
from itertools import count, 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 P(n, k = 0):
    if n < 0:
        return 0
    elif n == 0:
        return 1
    elif k == 0 or k > n:
        k = n
    
    if (n, k) in memo:
        return memo[(n,k)]
    
    s = sum(P(n - p, p) for p in takewhile(
                lambda p: p <= min(n, k), gen_primes()))
    memo[(n,k)] = s
    return s

N = 5000
memo = { }
primes = [ 2 ]
a = head(ifilter(lambda a: a[1] > N, ((n, P(n)) for n in count(0))))
print a[1]

# Code2
from itertools import imap, count, ifilter

def solve(max_n):
    def mul(a, p):
        return [ sum(imap(lambda l: a[k-l], xrange(0, k + 1, p)))
                                        for k in xrange(max_n + 1) ]
    
    def is_prime(n):
        return all(n % p != 0 for p in xrange(2, n))
    primes = filter(is_prime, xrange(2, max_n + 1))
    
    a = [ 1 ] + [ 0 ] * max_n
    for p in primes:
        a = mul(a, p)
        if a[p] > N:
            return p
    
    return 0

N = 5000
for n in ifilter(lambda n: n > 0, imap(solve, (1 << e for e in count(1)))):
    print n
    break

# Code3
from itertools import ifilter, imap, takewhile, count, islice
import sys

def gen(x):
    a, next, args = x
    for k in count():
        if k < len(a):
            yield a[k]
        else:
            e = next(a, args)
            a.append(e)
            yield e

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

def last(a):
    for e in a:
        pass
    return e

sys.setrecursionlimit(10000)
primes = [ 2 ]
g_prime = (primes,
        lambda a, d: head(ifilter(lambda n: all(imap(lambda p: n % p != 0,
        takewhile(lambda p: p * p <= n, primes))), count(a[-1] + 1))), -1)

def mul(g, p):
    def coeff(b, n):
        if n < 0:
            return 0
        
        x = b.next()
        if n % p == 0:
            return x + coeff(b, n - 1)
        else:
            return coeff(b, n - 1)
    
    def next(a, args):
        return coeff(gen(g), len(a))
    
    dic_ary[p] = [ ]
    return (dic_ary[p], next, p)

dic_ary = { }
dic_ary[0] = [ 1 ]
g = (dic_ary[0], lambda n, a: 0, 0)
for p in gen(g_prime):
    g = mul(g, p)
    if last(islice(gen(g), p + 1)) - 1 > 5000:
        print p
        break