Project Euler 123

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

(p_n-1)ⁿ + (p_n+1)ⁿ = 2n p_n (n : 奇数)
(p_n-1)ⁿ + (p_n+1)ⁿ = 2 (n : 偶数)

なので、あとは素数を小さいものから順に出すだけです。上限が決まっていないので、エラトステネスのふるいを少しずつしていくのが速いです。

from itertools import *

def gen_primes():
    a = [ True ] * L
    for p in (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 p in primes:
        yield p
    
    for m in count(1):
        start = m * L
        end = start + L
        a = [ True ] * L
        for p in takewhile(lambda p: p * p < end, primes):
            k0 = (start + p - 1) / p * p - start
            for k in xrange(k0, L, p):
                a[k] = False
        
        for p in (start + k for k in xrange(L) if a[k]):
            primes.append(p)
            yield p

def mod(n, p):
    return 2 if n % 2 == 0 else 2 * n * p % (p * p)

N = 10 ** 10
L = 10000
print next(n for n, p in izip(count(1), gen_primes()) if mod(n, p) > N)