Project Euler - Problem 196

昨日はみんな大好き HUB の日 ということで鳥貴族に行ってきました．それはそれとして Problem 196 をやります．

Prime triplets

Problem 196 (日本語訳)



だるいので和訳やめた．日本語訳のリンクを貼ったのでそのようにしてください．

自然数を三角形状に並べたとき，n 行目に存在する素数のうち prime triplet の要素であるもの総和を S(n) とする．この時の S(5678027) + S(7208785) を求める問題．

素数を求める段階と prime triplet を判定する段階の 2 ステップあります．



1. 素数を求める

素数の列挙といえばエラトステネスの篩ですが，7208785 行目の右端は 25983294192505 > 10¹³ なのでナイーブに実装するとメモリが絶滅します．後述しますが実のところ必要なのは n 行目を中心とした 5 行の情報だけなので，区間を限定することで空間計算量を削減する区間ふるいという方法が使えます．

//区間ふるい : [L, R) の素数を求める
bool* sieve(lint L, lint R) {
    

    // √R までの素数を求めておく
    lint sqn = (lint)sqrt(R) + 1;
    bool* pre = new bool[sqn];
    for ( lint i = 0; i < sqn; i++ ) {
        pre[i] = true;
    }
    pre[0] = pre[1] = false;
    for ( lint i = 2; i * i < sqn; i++ ) {
        if ( pre[i] ) {
            for ( lint j = i * i; j < sqn; j += i ) {
                pre[j] = false;
            }
        }
    }

    lint width = R - L;

    // 後の 8 近傍探索のために余白をとっておく
    lint l_margin = 2;
    lint r_margin = 5;
    bool* p = new bool[ l_margin + width + r_margin ];
    for ( lint i = 0; i < l_margin; i++ ) p[i] = false;
    for ( lint i = l_margin + width; i < l_margin + width + r_margin; i++ ) p[i] = false;
    p += l_margin;

    // √R までの素数の倍数をふるい落とす
    for ( lint i = 0; i < width; i++ ) p[i] = true;
    for ( lint i = 0; i < sqn; i++ ) {
        if ( pre[i] ) {
            lint q;
            if ( L % i == 0 ) q = L / i;
            else q = L / i + 1;
            for ( lint j = i * q; j < R; j += i ) {
                p[j - L] = false;
            }
        }
    }
    
    delete[] pre;
    
    return p;
}

2. prime triplet を判定する

y 行目 x 列の数 P = y(y - 1)/2 + x + 1 が素数であるとき，それが prime triplet の要素であるため条件は，

• P を含む素数の 8 連結 な連結領域の要素が 3 以上である

と言い換えることができる．これを基に近傍探索を行う． 18 ~ 24 行の様子です．

素数 P の 8 近傍について，

• 素数が 2 個以上 : prime triplet の要素 ( P = 199, 211 など )
• 素数が 0 個 : prime triplet の要素ではない ( P = 223, 239 など )
• 素数が 1 個の場合(これを Q とおく)，Q の 8 近傍について，
  • 素数が 2 個以上 : prime triplet の要素 ( (P, Q) = (197, 179), (233, 211) など )
  • それ以外 : prime triplet の要素ではない ( (P, Q) = (229, 251), (241, 263) など )

これによって prime triplet の要素かどうか判定できる．素数の配列は 5 行だけ確保しておけばいいことがここからわかる．

// prime triplet 判定
bool is_prime_triplet(bool** p, lint x, lint y) {

    // 3 2 1
    // 4   0
    // 5 6 7
    static const lint dx[8] = { 1,  1,  0, -1, -1, -1,  0,  1 };
    static const lint dy[8] = { 0, -1, -1, -1,  0,  1,  1,  1 };

    lint count = 0;
    lint xtmp, ytmp;
    
    for ( lint dir = 0; dir < 8; dir++ ) {
        if ( p[ y + dy[dir] ][ x + dx[dir] ] ) {
            count++;
            if ( count == 1 ) {
                xtmp = x + dx[dir];
                ytmp = y + dy[dir];
            }
        }
    }

    if ( count >= 2 ) return true;
    if ( count == 0 ) return false;

    count = 0;

    for ( lint dir = 0; dir < 8; dir++ ) {
        if ( p[ ytmp + dy[dir] ][ xtmp + dx[dir] ] ) {
            count++;
        }
    }

    if ( count >= 2 ) return true;

    return false;
}

変位を配列に格納して for で回すやつはよくあるアレですね．



以下ソースコード全体

#include <iostream>
#include <cmath>
#include <time.h>

using namespace std;
typedef long long int lint;

//区間ふるい : [L, R) の素数を求める
bool* sieve(lint L, lint R) {
    

    // √R までの素数を求めておく
    lint sqn = (lint)sqrt(R) + 1;
    bool* pre = new bool[sqn];
    for ( lint i = 0; i < sqn; i++ ) {
        pre[i] = true;
    }
    pre[0] = pre[1] = false;
    for ( lint i = 2; i * i < sqn; i++ ) {
        if ( pre[i] ) {
            for ( lint j = i * i; j < sqn; j += i ) {
                pre[j] = false;
            }
        }
    }

    lint width = R - L;

    // 後の 8 近傍探索のために余白をとっておく
    lint l_margin = 2;
    lint r_margin = 5;
    bool* p = new bool[ l_margin + width + r_margin ];
    for ( lint i = 0; i < l_margin; i++ ) p[i] = false;
    for ( lint i = l_margin + width; i < l_margin + width + r_margin; i++ ) p[i] = false;
    p += l_margin;

    // √R までの素数の倍数をふるい落とす
    for ( lint i = 0; i < width; i++ ) p[i] = true;
    for ( lint i = 0; i < sqn; i++ ) {
        if ( pre[i] ) {
            lint q;
            if ( L % i == 0 ) q = L / i;
            else q = L / i + 1;
            for ( lint j = i * q; j < R; j += i ) {
                p[j - L] = false;
            }
        }
    }
    
    delete[] pre;
    
    return p;
}

//n 段目の素数判定
bool* triangle_sieve(lint n) {
    return sieve(n * (n - 1) / 2 + 1, n * (n + 1) / 2 + 1);
}

//n 段目を中心とした高さ 2 * h + 1 段の素数判定
bool** multi_triangle_sieve(lint n, lint h) {
    bool** p = new bool*[2 * h + 1];
    p += h;
    for ( lint i = -h; i <= h; i++ ) {
        p[i] = triangle_sieve(n + i);
    }
    return p;
}

// prime triplet 判定
bool is_prime_triplet(bool** p, lint x, lint y) {

    // 3 2 1
    // 4   0
    // 5 6 7
    static const lint dx[8] = { 1,  1,  0, -1, -1, -1,  0,  1 };
    static const lint dy[8] = { 0, -1, -1, -1,  0,  1,  1,  1 };

    lint count = 0;
    lint xtmp, ytmp;
    
    for ( lint dir = 0; dir < 8; dir++ ) {
        if ( p[ y + dy[dir] ][ x + dx[dir] ] ) {
            count++;
            if ( count == 1 ) {
                xtmp = x + dx[dir];
                ytmp = y + dy[dir];
            }
        }
    }

    if ( count >= 2 ) return true;
    if ( count == 0 ) return false;

    count = 0;

    for ( lint dir = 0; dir < 8; dir++ ) {
        if ( p[ ytmp + dy[dir] ][ xtmp + dx[dir] ] ) {
            count++;
        }
    }

    if ( count >= 2 ) return true;

    return false;
}

lint S(lint n) {

    lint ret = 0;
    lint begin = n * (n - 1) / 2 + 1;

    lint h = 2;
    bool** p = multi_triangle_sieve(n, h);

    for ( lint x = 0; x < n; x++ ) {
        if ( p[0][x] ) {
            if ( is_prime_triplet(p, x, 0) ) {
                ret += begin + x;
            }
        }
    }

    return ret;
}

void PE0196()
{
    cout << S(5678027) + S(7208785) << endl;
}

int main() {

    clock_t start, end;
    start = clock();

    PE0196();

    end = clock();
    printf("%d msec.\n", end - start);

    return 0;
}

なんとなく C++ のほうが速いかな？と思ったのでそのようにしました．4 秒くらい．