概素数の列挙(1)

$k$-概素数の列挙、それなりに頑張って高速化したので手法について記録を残しておきます。

概素数とは

ある自然数$k$に対して自然数$n$が$k$-概素数($k$-almost prime)であるとは次のように定義されます。

$k$個の素数${p_1, p_2, \dots, p_k}$(それぞれの素数は互いに等しくてもよい)が存在して、$n = p_1 p_2 \dots p_k$が成り立つこと。^{1 2}
全く同じことですが、$\Omega (n)$をnの素因数の個数と定義したときに$\Omega (n) = k$であること、と表すこともできます。

例えば、$k=3$であるときの$k$-概素数を小さいものから順に並べていくと
${8 = 2^3, \quad 12 =2^2 \cdot 3,\quad 18 = 2 \cdot 3^2, \quad 20 = 2^2 \cdot 5, \quad 27 = 3^3, \quad \dots}$
となります。定義からわかる通り、もし$k=1$であれば$k$-概素数は素数と、もし$k=2$であれば$k$-概素数は半素数と一致します。また、この記事では$k=0$に対する$k$-概素数は$1$ただ一つのみであると定義することにします。

目標

この記事においての目標は、任意に与えられた自然数$k, n$に対して$n$未満の$k$-概素数を昇順の配列で返すこと、その処理を出来る限り高速化することです。C++20と、コンパイラにはgccを使用しコンパイラオプションで-O2の最適化をおこないます。今回は並列処理は明示的には行わない方針でいきます。

記号の定義

これからの議論を簡単にするためいろいろ数学記号を定義していきます。とは言っても実際に使用するのは大抵は初めに紹介する一つだけだと思うので、それ以外は無視してしまってもあまり問題はありません。

$\Omega (n)$
　nの素因数の個数。

$m \mid n$
　$m$が$n$を割り切る。

$m \nmid n$
　$m$が$n$を割り切らない。

$P_k$
　$k$-概素数の集合。

$P_k^n$
　$n$未満の$k$-概素数の集合。

アルゴリズム Ver 1

今回は、まずは最も簡単なアルゴリズムから紹介。
任意の自然数$n$と任意の非負整数$i$について、$0 \leqq i < \Omega(n)$であることと、次の三条件を満たす自然数$m$が存在することは同値となります。

1. $m < n$
2. $m$は$n$を割り切る
3. $\Omega(m) = i$

自然数$n$に対して、すでに$n$未満の自然数が全て素因数の個数によって$P[0], \ P[1], \ P[2], \ \dots \ $に分類されているという状況を考えてみましょう($P[k]$が$k$-概素数の集合に対応します)。そこから$n$の素因数の個数を求める方法は上に述べた事実から次のようなものが考えられます。

1. $i = 0$をとる
2. $P[i]$の中に$n$を割り切るような自然数があるかを探す
3. もし条件を満たす自然数が存在すれば$i$を1増やして 2. に戻る。条件を満たす自然数が存在しなければ$n$を$P[i]$に格納する。

ただし、今回の目標は$k$-概素数を求めることなので、もし$i$が$k+1$以上になった場合の$n$はどこにも格納せず捨ててしまっても構いません。

これをプログラムに落とし込むと以下のようになります。

std::vector<int64_t> almprm1(const int k, const int64_t n) {

    std::vector<std::vector<int64_t>> P(k+1);
    std::vector<int64_t> Pk;

    for (int64_t i = 1; i < n; ++i) {
        for (int64_t j = 0; j <= k; ++j) {
            if (std::ranges::any_of(P[j], [i](const int64_t x){return i%x == 0;})) {
                continue;
            }
            P[j].push_back(i);
            break;
        }
    }

    /* 二次元配列の一部にRVOが働くのかよくわからないのでとりあえずmoveを挟む */
    Pk = std::move(P[k]);

    return Pk;
}

さて、このアルゴリズム、実は実装が滅茶苦茶簡単なだけで速度はかなり遅いです。自然数$i$に対して割り切れるか確認するのは$i/2$までで良いとか色々改善の余地はあるんですが、それらを全部組み合わせてもどうしようもなく遅いです。ということで発想を転換しましょう。このアルゴリズムでは$n$未満の自然数を素因数の個数で振り分けるということをしましたが、そうではなく添え字$i$番目に$i$の素因数の個数を格納するような配列$PF$を作成し、そこから$k$-概素数の集合を取り出すという方針で別のアルゴリズムを作ります。³

アルゴリズム Ver 2.1

新しいアルゴリズムは以下のようなものです。

1. 長さ$n$の配列$PF$を用意し0番目は-1、それ以外の要素は全て0で初期化する。
2. $i = 2, \dots, n-1$まで3. を繰り返す。そのあと5. へ。
3. $PF[i] \neq 0$であれば何もしない、$PF[i] = 0$であれば$r = i, i^2, i^3, \dots, i^j$($j$は$i^j < n, \quad j \leqq k+1$を両方満たすような最大の整数)それぞれに対して4. を行う。
4. 配列$PF$の${ \ r, \ 2r , \ 3r, \ \dots \ }$番目の値を1増やす。
5. $PF[j] = k$となるような$j$の配列を返す。

プログラムにしたものがこちら。

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

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

    for (int64_t i = 2; i < n; ++i) {
        if (PF[i]) {
            continue;
        }
        for (int64_t j = 1, r = i; r < n && j <= k+1; ++j, r *= i) {
            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;
}

次回はこのアルゴリズムをさらに高速化していきます。

計算時間

単位は秒です。Ver 1とVer 2.1でn, kの範囲が異なることに注意。

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

1. $n$が$K$-概素数であることを$\Omega(n) \leqq K$であると定義することもあります。 ↩

2. 概素数のWikipedia ↩

3. 少しふわっとした話になってしまいますが、ここでの古いアルゴリズムと新しいアルゴリズムの関係は、素数列挙における試し割り法とエラトステネスの篩の関係によく似ています。 ↩