Project Euler - Problem 291

Panaitopol Primes

素数 $p$ がある正の整数 $x, y$ に対して $\displaystyle p = \frac{ x^{4} - y^{4} }{ x^{3} + y^{3} } $ を満たすとき，$p$ を Panaitopol 素数と呼ぶ．

$5×10^{15}$ 未満の Panaitopol 素数はいくつあるか．

(日本語訳より)

Difficulty rate 45% は詐欺

小さい値で実験をする

f:id:my316g:20180206035802p:plain

結果より，

$d = \mathrm{GCD}(x, y)$
$x = d x_{0}$
$y = d y_{0}$

として，

$x - y = d$
$x_{0} - y_{0} = 1$
$p = 2 d - 1$

であるらしいことがわかる．（証明は知らん）

で，これより，

$\displaystyle p = \frac{ x^{4} - y^{4} }{ x^{3} + y^{3} } = \cdots = \frac{ d (x_{0} - y_{0})(x_{0}^{2} + y_{0}^{2}) }{ x_{0}^{2} - x_{0} y_{0} + y_{0}^{2} }$

$\displaystyle \quad = \frac{ d \cdot 1 \cdot \{ (x_{0} - y_{0})^{2} + 2 x_{0} y_{0} \} }{ (x_{0} - y_{0})^{2} + x_{0} y_{0} } = \frac{ d \cdot 1 \cdot ( 1^{2} + 2 x_{0} y_{0} ) }{ 1^{2} + x_{0} y_{0} }$

$\displaystyle \quad = \frac{ d(1 + 2 x_{0} y_{0}) }{ 1 + x_{0} y_{0} } = 2 d - 1 $

最後の等号は $p = 2d - 1$ を使った．整理すると，

$ d = 1 + x_{0} y_{0} $

という関係が得られ，

$ p = 2d - 1 = 2(1 + x_{0} y_{0}) - 1 = 2 x_{0} y_{0} + 1 $
$ \quad = 2 y_{0} (y_{0} + 1) + 1 \ \ (∵ \ x_{0} - y_{0} = 1 )$

$ \quad =y_{0}^{2} + (y_{0} + 1)^{2} $

がわかる．

あとは $y_{0}$ をインクリメントして $p$ が素数かつ $5 \times 10^{15}$ 以下であるものをカウントしていけばいい．

というのは嘘で，$N = 5 \times 10^{15}$ とすると $y_{0}$ は $\sqrt{N/2}$ 個位あって，素数判定 (試し割り) に $O(\sqrt{N})$ 掛かるので全体で $O(N)$ となり 1 分ルールを考えると余裕で TLE．(Difficulty rate が低いのは TLE 解でゴリ押ししている人間が一定数いるためだと思う)

あれこれ考えた結果，$2 n^{2} + 2 n + 1$ 型の素数をエラトステネスの篩っぽく求める方法をここなどを参考にして実装した．

直前の投稿に纏めておいたので興味のある人は見てください．

chillbrains.hateblo.jp

#include "bits/stdc++.h"
using namespace std;
typedef long long ll;

vector<ll> quadraticPrimes(ll N)
{
    vector<ll> v;
    vector<ll> ps;

    ll p;
    for (ll n = 0; (p = 2 * n * n + 2 * n + 1) <= N; n++) v.push_back(p);

    for (int i = 1; i < v.size(); i++) {
        p = v[i];
        if (p != 1) {
            if (i == 1 || p > ps.back()) ps.push_back(p);
            for (int j = i + p; j < v.size(); j += p) {
                while (v[j] % p == 0) v[j] /= p;
            }
            for (int j = p - (i + 1); j < v.size(); j += p) {
                while (v[j] % p == 0) v[j] /= p;
            }
        }
    }

    return ps;
}

int main()
{
    clock_t start, end;
    start = clock();

    //cin.tie(0);
    //ios::sync_with_stdio(false);

    const ll N = 5000000000000000LL;

    cout << quadraticPrimes(N).size() << endl;

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

    return 0;
}

10.7 sec くらい．

絶滅

どうでもいい

Project Euler - Problem 291

Panaitopol Primes