概素数の列挙(2)

前回の記事の続き。気が付いたら2か月経っていた。今回は前回の最後に紹介したアルゴリズムを、$k$-概素数の$k$が大きい場合に焦点を当ててさらに高速化していきます。

アルゴリズム Ver 2.2

前回最後に紹介したアルゴリズムでは初めから直接的に素因数の個数の配列$PF$を作成していましたが、今回は素数の配列$P$をエラトステネスの篩で作成しそれを基に素因数の個数の配列$PF$を作成するようにアルゴリズムを変更します。この変更自体は処理の高速化にはなりません。これから先の高速化をしていくための過渡的なものになります。
エラトステネスの篩それ自体の実装方法ですが、大変なのでここでは解説しません。その代わりですがこちらの記事を参考にするとよいでしょう。この記事のエラトステネスの篩実装は自分の知る限りではQiita上最速です。ただしかなり内容は難しいので注意。この記事中で使用するエラトステネスの篩の実装はこちらの記事の(1), (2), (6)の高速化を適用させたものを使用していきます。

std::vector<int64_t> almprm2_2(const int8_t k, const int64_t n) {

    /* n未満の素数の配列 */
    const std::vector<int64_t> P = sieve(n);

    std::vector<int8_t> PF(n);
    PF[0] = -1;

    for (auto p : P) {
        for (int64_t j = 1, r = p; r < n && j <= k+1; ++j, r *= p) {
            for (int64_t l = r; l < n; l += r) {
                ++PF[l];
            }
        }
    }

    std::vector<int64_t> Pk(std::ranges::count(PF, k));
    auto Pk_iter = std::begin(Pk);
    for (int64_t j = 0; j < n; ++j) {
        if (PF[j] == k) {
            *(Pk_iter++) = j;
        }
    }

    return Pk;
}

さて、このソースコードに対してこれから本命の高速化を行っていきます。
このソースコードでははじめに$n$未満の素数配列$P$を作成していますが、この部分には無駄があります。例として$k=3, n=100$の場合を考えてみましょう。
この場合上のソースコードだと素数配列$P$は${2, \ 3, \ 5, \ \dots \ 23, \ 27, \dots 83, \ 89, \ 97}$となりますが、実はこれらの値たちのうち27以降は全て必要ありません。27以上の素数に対してそれらを素因数に含む$k$-概素数はすべて100以上であり最終的な出力結果には含まれないからです(ex. $27$を素因数に含む$k$-概素数の最小は$2^2 \cdot 27 = 108$)。一般化した$k, n$について考えてみると$P$は$\left\lceil \frac{n}{2^{k-1}} \right\rceil$未満の素数配列とすれば十分となります（ただし、$\lceil x \rceil$の定義はx以上の整数のうちで最小のもの）（ただし$k = 0$の場合バグるので条件分岐をしましょう）。これだけで随分と高速になります。

もう一つ、素因数の個数の配列$PF$の作成の部分にも無駄があり、高速化することができます。
$k$-概素数のうち最少の値がいくつなのかを考えると、それは$2^k$だと容易に分かります。これは逆に言うと$2^k$未満の整数は必ず$k$-概素数にならないということなので、$PF$の$2^k$未満に対応する部分は始めから飛ばし、$PF[i]$がもともと$i$の素因数の個数を表していたものから$i+2^k$の素因数の個数を表すようアルゴリズムを変換することでさらに処理を高速化できます。
そのためには$PF$の確保するサイズを$2^k$分減らし（$n$が$2^k$未満である場合は場合分けをしましょう）、ソースコード中for文の変数$l$の初期値を、$r$の倍数で$2^k$以上になるような最も小さい値から$2^k$を引いた値である$r - (2^k - 1) \% r - 1$に、あとはその他の部分も適切にfor分の範囲を変更しておきましょう。
各$k$に対して$PF$のサイズがいくら減らせるかを考えると明らかですが、この高速化は$k$が小さいときにはほぼ効果はありません。むしろ計算が複雑化する分遅くなることもあるのでそこだけ注意。

そうして出来上がったソースコードがこちらです。

std::vector<int64_t> almprm2_2(const int8_t k, const int64_t n) {

    const int64_t kp = 1 << k;
    if (k == 0) [[unlikely]] return {1};
    if (n <= kp) [[unlikely]] return {1};

    const auto P = sieve(-(-n >> (k-1)));

    std::vector<int8_t> PF(n-kp);

    for (auto p : P) {
        for (int64_t j = 1, r = p; r < n && j <= k+1; ++j, r *= p) {
            for (int64_t l = r - (kp-1)%r - 1; l < n-kp; l += r) {
                ++PF[l];
            }
        }
    }

    std::vector<int64_t> Pk(std::ranges::count(PF, k));
    auto Pk_iter = std::begin(Pk);
    for (int64_t j = 0; j < n-kp; ++j) {
        if (PF[j] == k) {
            *(Pk_iter++) = j + kp;
        }
    }

    return Pk;
}

次回は全く別のアルゴリズムを紹介します。

計算時間

比較のため前回のものも記載、カッコで記載しているのは計算時間全体のうち、エラトステネスの篩にかけた時間。
単位は秒。

Ver 1

	k = 1	k = 2	k = 4	k = 8
n = 10^3	-	0.002	0.003	0.003
n = 10^4	0.027	0.149	0.309	0.332
n = 10^5	1.437	6.881	14.569	22.461

Ver 2.1

	k = 1	k = 2	k = 4	k = 8	k = 16
n = 10^5	0.002	0.002	0.002	0.002	0.002
n = 10^6	0.026	0.033	0.033	0.025	0.026
n = 10^7	0.327	0.404	0.403	0.330	0.328
n = 10^8	4.687	5.185	5.492	4.890	4.584

（この記事１つ目のソース、比較用）

	k = 1	k = 2	k = 4	k = 8	k = 16
n = 10^5	0.003(-)	0.003(-)	0.003(-)	0.002(-)	0.002(-)
n = 10^6	0.027(0.002)	0.034(0.002)	0.033(0.002)	0.026(0.002)	0.024(0.002)
n = 10^7	0.322(0.028)	0.386(0.028)	0.383(0.027)	0.320(0.029)	0.295(0.027)
n = 10^8	4.722(0.297)	5.362(0.343)	5.450(0.332)	4.825(0.333)	4.698(0.290)

Ver 2.2¹

	k = 1	k = 2	k = 4	k = 8	k = 16
n = 10^5	0.003(-)	0.003(-)	0.003(-)	0.002(-)	0.001(-)
n = 10^6	0.032(0.002)	0.032(0.001)	0.029(-)	0.021(-)	0.012(-)
n = 10^7	0.348(0.028)	0.373(0.013)	0.303(0.002)	0.226(-)	0.188(-)
n = 10^8	5.018(0.311)	5.275(0.160)	5.185(0.036)	3.890(0.001)	2.919(-)

補足

$\lceil \frac{n}{2^{k-1}} \rceil$は以下のように求められます。

-(-n >> (k-1))

ただしこれはnがsignedである場合にのみ成り立つ式なので注意。

1. この表の通り、実用的には$k=1$の場合にエラトステネスの篩の値をそのまま返すようにすれば大幅な速度アップになります。今回は面倒だったのであくまでもアルゴリズムの部分に焦点を当てているのでそういった高速化はしませんでした。 ↩