Project Euler 133

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

素数pでR(10ⁿ)が割り切れるということは、A(p)でR(10ⁿ)が割り切れるということです。そのようなnが存在するということは、A(p)は2か5以外の因子を含まないということです。

from itertools import *
from fractions import gcd

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

def factorize(max_n):
    a = range(max_n + 1)
    b = [ () for n in a ]
    for p in takewhile(lambda p: p * p <= max_n,
                ifilter(lambda n: a[n] == n, xrange(2, max_n + 1))):
        for n in xrange(p * 2, max_n + 1, p):
            e, a[n] = div_pow(a[n], p)
            b[n] = b[n] + ((p, e),)
    
    for n in ifilter(lambda n: a[n] > 1, xrange(2, max_n + 1)):
        b[n] = b[n] + ((a[n], 1),)
    
    return b

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

def gen_divisors(f, k = 0):
    if k == len(f):
        yield 1
    else:
        p, e0 = f[k]
        for n in gen_divisors(f, k + 1):
            yield n
            for e in range(1, e0):
                yield p ** e * n
            yield p ** e0 * n

def get_exp(n, p):
    return 0 if n % p != 0 else 1 + get_exp(n / p, p)

def A(n):
    f = fs[n-1]
    for d in sorted(list(gen_divisors(f))):
        if pow(10, d, n) == 1:
            e3 = get_exp(n, 3) - get_exp(d, 3)
            if e3 > 0:
                return d * 3 ** e3
            else:
                return d

def is_prime(n):
    f = fs[n]
    return len(f) == 1 and f[0][1] == 1

def is_valid(n):
    def g(f):
        if len(f) == 0:
            return True
        
        p = f[0][0]
        if p == 2 or p == 5:
            return g(f[1:])
        else:
            return False
    
    if gcd(n, 30) != 1:
        return False
    return g(fs[A(n)])

N = 10 ** 5
fs = factorize(N - 1)
print sum(n for n in xrange(2, N) if is_prime(n) and not is_valid(n))