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
e1 ≥ e2 ≥ … ≥ 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