n^2+1型の素数の列挙

n² + 1という形の素数を求める問題はProject Eulerで頻出なので簡単なアルゴリズムを解説します。

例として200以下のこの形の素数を列挙します。
まず、n² + 1を順に並べます。

2 5 10 17 26 37 50 65 82 101 122 145 170 197

2で割り切れる数が1つおきにあるので、それを最初だけ残して2で割ります。

2 5 5 17 13 37 25 65 41 101 61 145 85 197

右に移動して5です。1でないので素数です（後述）。これをpとします。

n = 2
n² + 1 ≡ 0 (mod p)

そうすると、これもpで割り切れることが分かります。

(n + kp)² + 1

また、これもpで割り切れます。

(kp - n)² + 1

すなわち、n = 7, 12, 17, ..., 3, 8, 13, ...も5で割り切れます。これらを5で割れるだけ割ります。

2 5 1 17 13 37 1 13 41 101 61 29 17 197

次は1なのでその右の17について同様に処理します。n = 4なので、17 - n = 13について処理します。

2 5 1 17 13 37 1 13 41 101 61 29 1 197

これを最後まで繰り返します。

2 5 1 17 13 37 1 1 41 101 61 29 1 197

一度も割られていないものがn²+1型の素数となります。

2 5 17 37 101 197

from itertools import ifilter
from math import sqrt

def div_pow(n, d):
    return n if n % d != 0 else div_pow(n / d, d)

def sieve(N):
    M = int(sqrt(N - 1))
    a = [ n * n + 1 for n in xrange(M + 1) ]
    for n in xrange(3, M + 1, 2):
        a[n] /= 2
    
    for p, n in ((a[n], n) for n in xrange(2, M + 1) if a[n] != 1):
        for q in xrange(n + p, M + 1, p):
            a[q] = div_pow(a[q], p)
        for q in xrange(p - n, M + 1, p):
            a[q] = div_pow(a[q], p)
    
    return [ a[n] for n in xrange(1, M + 1) if a[n] == n * n + 1 ]

print sieve(200)

1でなければ素数というのは、もし別に素因数があったとするとその前に出てきているはずからです。nに対して素因数pがあって、

n² + 1 ≡ 0 (mod p)
p < n

となります（複数の素因数のうちのどれかがこうなります）。そうすると、(n - p)² + 1もpで割り切れなくてはならなくなります。

Project Eulerに必要な用語集
http://d.hatena.ne.jp/inamori/20091225/p1