Project Euler 146

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

すぐにn²は30で割って余りが10でなければならないとわかるでしょう。すなわち、nは30で割って余りが10か20です。しかし、このnに対してまともに素数判定をやっていてはいくら時間があっても足りません。
そこでフェルマーテストを行います。これで素数であることはわかりませんが、合成数であることはわかるので、かなり振り落とせるはずです。実際にほぼ正しい答えしか出てこないようです。
しかし、このほんのわずかになった解の候補にさえ、ちゃんと素数判定しようとすると時間がかかります。最初は6で割って1か5余る数で割っていけばいいだろうと思っていたのですが、時間がかかるのでエラトステネスのふるいで必要な素数を求め、さらにn² 〜 n² + 27でエラトステネスのふるいをかけたのですが、それほど速くなりませんでした。結局答えが出るまでに20分ほどかかりました。答えを出すだけならミラー・ラビンでもいいのでしょうが。

from itertools import *

c = [ 1, 3, 7, 9, 13, 27 ]
d = [ 19, 21 ]

def sieve(max_n):
    a = [ True ] * L
    for p in takewhile(lambda n: n * n < L,
                    (n for n in xrange(2, L) if a[n])):
        for k in xrange(p * 2, L, p):
            a[k] = False
    primes = [ n for n in xrange(2, L) if a[n] ]
    
    for offs in xrange(L, max_n + 1, L):
        a = [ True ] * L
        for p in takewhile(lambda n: n * n < offs + L, primes):
            for m in xrange((offs + p - 1) / p * p, offs + L, p):
                a[m-offs] = False
        primes.extend(k + offs for k in xrange(L) if a[k])
    return primes

def sieve2(n):
    a = [ True ] * 28
    for p in takewhile(lambda p: p * p <= n + 27, primes):
        k0 = (n + p - 1) / p * p - n
        for k in xrange(k0, 28, p):
            a[k] = False
    return all(a[k] for k in c) and all(not a[k] for k in d)

def is_prime(n):
    return all(n % p != 0 for p in takewhile(lambda p: p * p <= n, primes))

def gen_n(N):
    return chain(xrange(10, N, 30), xrange(20, N, 30))

def Fermat_test(n):
    return pow(2, n - 1, n) == 1

def is_consecutive_primes(n):
    if n % 7 == 0 or n % 13 == 0:
        return False
    return all(imap(Fermat_test, ((n * n + e) for e in c))) \
                and sieve2(n * n)

N = 150 * 10 ** 6
L = 10 ** 5 + 1
primes = sieve(N + 4)
print sum(ifilter(is_consecutive_primes, gen_n(N)))