Project Euler 124

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

10万より小さい自然数を素因数分解し、radを計算してソートするだけです。
もう少し速い方法はないでしょうか。まず、べき乗を含まない自然数を出します。10までなら2,3,5,6,7,10です。例えばrad(n)が6になるnは、2^e3^fの形です。この10万より小さい自然数の個数を数えます。小さいほうから足していって1万以上になるまで繰り返します。そうなったら今のnについてもう一度rad(n)になる自然数を、今度は個数ではなく数を出します。それをソートしてちょうど1万になる数を得ます。これで最初の方法より1桁くらい速いです。

from itertools import *

def sieve(max_n):
    a = [ True ] * (max_n + 1)
    for p in (n for n in xrange(2, max_n + 1) if a[n]):
        for k in xrange(p * 2, max_n + 1, p):
            a[k] = False
    return [ n for n in xrange(2, max_n + 1) if a[n] ]

def gen_factors():
    def div_pow(n, p):
        e = 0
        while n % p == 0:
            e += 1
            n /= p
        return e, n
    
    start = 2
    while True:
        a = range(start, start + L)
        b = [ [] for k in xrange(start, start + L) ]
        for p in takewhile(lambda p: p * p < start + L, primes):
            if p >= start:
                begin = p * 2
            else:
                begin = (start + p - 1) / p * p
            for n in xrange(begin, start + L, p):
                k = n - start
                e, a[k] = div_pow(a[k], p)
                b[k].append((p, e))
        
        for k in xrange(L):
            n = start + k
            if len(b[k]) == 0:
                b[k] = [ (n, 1) ]
                primes.append(n)
            elif a[k] > 1:
                b[k].append((a[k], 1))
        
        for f in b:
            yield f
        
        start += L

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

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

def is_equal_rad(f):
    return all(g[1] == 1 for g in f)

def gen_pows(p):
    m = 1
    while True:
        yield m
        m *= p

def num_rad(f, n):
    if len(f) == 0:
        return 1
    
    p = f[0][0]
    return sum(num_rad(f[1:], n / m) for m in
                takewhile(lambda m: m <= n, gen_pows(p)))

def gen_rad(f, n):
    if len(f) == 0:
        yield 1
    else:
        p = f[0][0]
        for m in takewhile(lambda m: m <= n, gen_pows(p)):
            for q in gen_rad(f[1:], n / m):
                yield m * q

N = 100000
M = 10000
L = 1000
primes = sieve(L)
counter = 1
for f in ifilter(is_equal_rad, gen_factors()):
    n = num_rad(f, N / f_value(f))
    counter += n
    if counter >= M:
        break
k = n - counter + M - 1
print sorted(list(gen_rad(f, N / f_value(f))))[k] * f_value(f)