Project Euler 110(1)

http://projecteuler.net/index.php?section=problems&id=110

n = p1e1…pmem

とすると、解の個数は、

((2e1 + 1)…(2em + 1) + 1) / 2

なので、指数のみで決まります。例えば、2e1・3e2と2e1・5e2では、解の個数は同じでも前者の方が小さくなります。また、23・34より24・33のほうが小さくなります。k番目の素数pkとすると、答えは

p1e1…pmem
e1e2 ≥ … ≥ em

の形でなければならないことがわかります。この形の自然数を勝手にヤング数と呼ぶことにしましょう(ヤング図形からの類推)。ヤング数を昇順に並べられれば、最初に解の個数が400万を超えた数が答えになります。
最小のヤング数は2です。2をリストに入れて、ここから素数を追加してヤング数を作っていきます。2からは22と2・3が生成されます。2はリストから取り出して、22と2・3のうち小さいのは22の方です。これをリストから取り出して、23と22・3が生成されるので、これをリストに入れます。このようにすればヤング数を並べられます。
これで10秒でした。もっと速くならないのでしょうか。

from itertools import *

def get_primes(k):
    if k == len(primes):
        for n in count(primes[-1] + 1 if k > 0 else 2):
            if all(n % m != 0 for m in xrange(2, n)):
                primes.append(n)
                break
    return primes[k]

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

def gen_Young():
    def f_tuple(f):
        return (f, f_value(f))
    
    def gen_next(f0):
        if len(f0) == 1 or f0[-2][1] > f0[-1][1]:
            p, e = f0[-1]
            yield f_tuple(f0[:-1] + [(p, e + 1)])
        yield f_tuple(f0 + [(get_primes(len(f0)), 1)])
    
    def index(a, y, begin = 0, end = -1):
        if end == -1:
            if len(a) == 0:
                return 0
            elif y[1] > a[-1]:
                return len(a)
            end = len(a)
        
        if begin == end - 1:
            return end
        else:
            mid = (begin + end) / 2
            if a[mid][1] < y[1]:
                return index(a, y, mid, end)
            else:
                return index(a, y, begin, mid)
    
    a = [ f_tuple([(get_primes(0), 1)]) ]
    while True:
        x = a.pop(0)
        yield x
        for y in gen_next(x[0]):
            k = index(a, y)
            a = a[:k] + [ y ] + a[k:]

def count_solutions(f):
    return (reduce(lambda x, y: x * (y[1] * 2 + 1), f, 1) + 1) / 2

N = 4 * 10 ** 6
primes = []
for x in gen_Young():
    if count_solutions(x[0]) > N:
        print x[1]
        break