$2n^{2} + 2n + 1$ 型素数の列挙

Project Euler 291 を解いていて 2n² + 2n + 1 型の素数を列挙する必要性が生じたので纏めておく．

ここに n² + 1 の場合の方法がありそれを参考にしている．

なお証明は知らん

素数を格納する配列 p を用意する．

2n² + 2n + 1 型の整数を配列 v に並べる (1-based，1 行につき 10 個並べてある)．

5, 13, 25, 41, 61, 85, 113, 145, 181, 221,
265, 313, 365, 421, 481, 545, 613, 685, 761, 841,
925, 1013, 1105, 1201, 1301, 1405, 1513, 1625, 1741, 1861,
1985, 2113, 2245, 2381, 2521, 2665, 2813, 2965, 3121, 3281,
3445, 3613, 3785, 3961, 4141, 4325, 4513, 4705, 4901, 5101,
5305, 5513, 5725, 5941, 6161, 6385, 6613, 6845, 7081, 7321,
7565, 7813, 8065, 8321, 8581, 8845, 9113, 9385, 9661, 9941, ...

先頭から見ていく．v[1] = 5 に着目．5 は 1 でないので，素数配列 p に格納する．

5 の倍数の出現位置を調べると n = 1 + 5k, 5k - 2 (k = 1, 2, ...) の場合に限定されることがわかる．これを割れなくなるまで片っ端から 5 で割っていく．

5, 13, 1, 41, 61, 17, 113, 29, 181, 221,
53, 313, 73, 421, 481, 109, 613, 137, 761, 841,
37, 1013, 221, 1201, 1301, 281, 1513, 13, 1741, 1861,
397, 2113, 449, 2381, 2521, 533, 2813, 593, 3121, 3281,
689, 3613, 757, 3961, 4141, 173, 4513, 941, 4901, 5101,
1061, 5513, 229, 5941, 6161, 1277, 6613, 1369, 7081, 7321,
1513, 7813, 1613, 8321, 8581, 1769, 9113, 1877, 9661, 9941, ...

割った．次に v[2] = 13 を見る．13 は 1 でないので，素数配列 p に格納する．

13 の倍数の出現位置を調べると n = 2 + 13k, 13k - 3 (k = 1, 2, ...) の場合に限定されることがわかる．これを割れなくなるまで片っ端から 13 で割っていく．

5, 13, 1, 41, 61, 17, 113, 29, 181, 17,
53, 313, 73, 421, 37, 109, 613, 137, 761, 841,
37, 1013, 17, 1201, 1301, 281, 1513, 1, 1741, 1861,
397, 2113, 449, 2381, 2521, 41, 2813, 593, 3121, 3281,
53, 3613, 757, 3961, 4141, 173, 4513, 941, 29, 5101,
1061, 5513, 229, 457, 6161, 1277, 6613, 1369, 7081, 7321,
1513, 601, 1613, 8321, 8581, 1769, 701, 1877, 9661, 9941, ...

割った．次に v[3] = 1 を見るが，1 の場合はこれを無視する．次の v[4] = 41 を見る．41 は 1 でないので，素数配列 p に格納する．

41 の倍数の出現位置を調べると n = 4 + 41k, 41k - 5 (k = 1, 2, ...) の場合に限定されることがわかる (表はだるいので省略)．これを割れなくなるまで片っ端から 41 で割っていく．

次の v[5] = 61 も 1 でないので素数配列 p に格納し，n = 5 + 61k, 61k - 6 の位置に現れる 61 の倍数を割る．

このへんで流石に法則が見えてくる．つまり v[n] = p であったとすると v[n + pk], v[pk - (n + 1)] は p の倍数になっている．

それ以外の数はどうやら p の倍数ではないらしい．証明は適当にやってください．

5, 13, 1, 41, 61, 17, 113, 29, 181, 17,
53, 313, 73, 421, 37, 109, 613, 137, 761, 841,
37, 1013, 17, 1201, 1301, 281, 1513, 1, 1741, 1861,
397, 2113, 449, 2381, 2521, 1, 2813, 593, 3121, 3281,
53, 3613, 757, 3961, 101, 173, 4513, 941, 29, 5101,
1061, 5513, 229, 457, 101, 1277, 6613, 1369, 7081, 7321,
1513, 601, 1613, 8321, 8581, 29, 701, 1877, 9661, 9941, ...

で，次の v[6] = 17 について，これは直前に格納された素数 61 よりも小さいので素数配列 p に格納しない．ただ，n = 6 + 17k, 17k - 7 の位置に現れる 17 の倍数を割る操作は行う．

次の v[7] = 113 に対しても素数配列格納および除算の操作を行う．

…

この操作を最後まで行ったとき，素数配列 p には 2n² + 2n + 1 型の素数が格納されている．

typedef long long ll;

// N 以下の 2n^2 + 2n + 1 型素数の配列を返す
vector<ll> quadraticPrimes(ll N)
{
    vector<ll> v;
    vector<ll> ps;

    // N 以下の 2n^2 + 2n + 1 型整数を配列 v に格納
    ll p;
    for (ll n = 0; (p = 2 * n * n + 2 * n + 1) <= N; n++) v.push_back(p);

    // 配列 v の最後まで操作を行う
    for (int i = 1; i < v.size(); i++) {
        p = v[i];
        // 1 でない場合
        if (p != 1) {
            // 最後に素数配列に格納した素数より大きければ格納 (最初の 5 は格納する)
            if (i == 1 || p > ps.back()) ps.push_back(p);
            // 位置 n = i + kp にある p の倍数を p で割れるだけ割る
            for (int j = i + p; j < v.size(); j += p) {
                while (v[j] % p == 0) v[j] /= p;
            }
            // 位置 n = kp - (i + 1) にある p の倍数を p で割れるだけ割る
            for (int j = p - (i + 1); j < v.size(); j += p) {
                while (v[j] % p == 0) v[j] /= p;
            }
        }
    }

    return ps;
}