Project Euler 130

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

前回の方法でも十分速く答えが出ますが、もっと速い方法があります。前回は小さい長さのrepunitから一つずつ順に割り切れるかどうか調べました。しかし実際には全てを調べる必要はありません。

まず9が続く数を考えましょう。これは10^k - 1と表せます。これがnで割り切れるとするとします。

10^k ≡ 1(mod n)

GCD(n, 10) = 1だからオイラーの定理より、

10^φ(n) ≡ 1(mod n)

よって、上を満たす最小のkはφ(n)の約数となります。約数を列挙して昇順にソートして下から上を満たすかを調べればよいのです。べき乗の計算はバイナリ法が使えます。Pythonではこれを実装しなくてもpowという関数が用意されています。
この問題は9ではなく1の連続なので、nが3の倍数のときは注意が必要です。例えば、3は9を割りますが、1は割りません。長さを3倍してA(3) = 3となります。

from itertools import *
from fractions import gcd

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 factorize(m):
    def div_pow(n, p):
        e = 0
        while n % p == 0:
            e += 1
            n /= p
        return e, n
    
    while len(facts) <= m:
        start = len(facts)
        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) ]
                if start > 0:
                    primes.append(n)
            elif a[k] > 1:
                b[k].append((a[k], 1))
        
        facts.extend(b)
    
    return facts[m]

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 value_f(f):
    return reduce(lambda x, y: x * y[0] ** y[1], f, 1)

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

def phi(n):
    f = factorize(n)
    return reduce(lambda x, y: x * (y[0] - 1) * y[0] ** (y[1] - 1), f, 1)

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

def is_composite(n):
    f = factorize(n)
    return len(f) > 1 or f[0][1] > 1

N = 25
L = 10000
primes = sieve(L)
facts = []
g =ifilter(lambda n: gcd(n, 10) == 1 and (n - 1) % A(n) == 0,
                            ifilter(is_composite, count(2)))
print sum(islice(g, 25))