min25筛

ある積性関数 f(x) が与えられ、f(p) が素数 p に関する多項式で、かつ f(p^k) が高速に計算できる場合、sum_i=1ⁿ f(i) を求める。

例：n = 10¹⁰ に対して sum_i=1ⁿ σ₀(i^k) を求める。

ステップ1：素数の積性関数の前処理和

1 から n までの素数の個数を求める。

f(n, m) を、2 から n までの数のうち、素数であるか、または最小の素因数が p₁, ..., p_m (p_i は i 番目の素数) ではない数の個数と定義する。

すると、以下の漸化式が成り立つ：

これは、f(n, m-1) から遷移し、最小の素因数がちょうど p_m である数をふるい落とすことを意味する。ふるい落とす数を p_m で割り、残った数の最小の素因数は p_m, p_m+1, ... となる。ここで「素数」部分を削除することに注意する。

注意：p_m 自体をふるい落としたくないが、f は 1 を含まないため、後ろの部分で 1 をふるい落とすことはないので、要件を満たしている。

より一般的な場合。

同様に、f(n, m) を、2 から n までの数のうち、素数であるか、または最小の素因数が p₁, ..., p_m ではない数の関数和と定義する。このとき、合成数も素数と同じように「偽の」関数和を持つと考える。すると、以下の式が得られる：

（関数名が重複したため、F(n) で積性関数を表す；psum は事前に計算された素数関数の前処理和）

多項式であるため、F(p) = p^k の形式に分解できる。最終的にそれらを組み合わせる。こうして、偽の関数は完全な積性関数となり、F(p_m) を直接取り出すことができる。

• 複雑度の証明

関連するブログ記事を参照。

p_m² ≤ n であることが分かる。

積分によって証明でき、複雑度は O(n^0.75 / log n) となる。

ステップ2：すべての数の積性関数の和

g(n, m) を、1 から n までの数のうち、最小の素因数が p_m より大きい数の（真の）積性関数の和と定義する。

すると、以下の式が成り立つ：

これは、素数と合成数の寄与を分けて考えることを意味する。素数の寄与はすでに求めた f(n, ∞) - psum(m) である。合成数の場合、最小の素因数とその指数を列挙する必要がある。つまり、k > m の p_k を列挙し、さらに指数 e を列挙する。積性関数であるため、F(p_k^e) を直接取り出すことができる。p^k (k > 1) のような数は、素数の部分で列挙されていないし、合成数の部分で p_k^e を取り出して列挙された後も列挙されていないことに注意が必要である。そのため、[e ≠ 1] を加える必要がある。

複雑度：証明できない。状態数は上記と同様に O(n^0.75 / log n) であるが、遷移の複雑度はやや難解で、総複雑度は O(n^0.75 / log n) とされることもあれば、O(n^{1 - ε}) とされることもある。

洲閣筛

機能は基本的に同じだが、min25筛ほど高速ではない。しかし、一度にすべての ⌊n/i⌋ の前処理和の値を求めることができる。

（通常の洲閣筛とは少し異なり、これは hs-black から教わったものである）

ステップ1：素数の積性関数の前処理和

min25筛と完全に一致するため、f(n, m) をそのまま使用する。

ステップ2：すべての数の積性関数の和

状態は少し異なり、ステップ1の f 配列に似ている。

g(n, m) を、2 から n までの数のうち、素数または最小の素因数が p_m より大きい数の真の積性関数の和と定義する。

すると、以下の式が成り立つ：

これは、p_m を含まない部分は g(n, m) から直接遷移し、p_m を含む部分は指数 c を列挙する。p_m^c の形であれば、答えは F(p_m^c) となる。そうでなければ、p_m より大きい他の素因数が存在するため、答えは F(p_m^c)(g(⌊n / p_m^c⌋, m) - psum(m)) となる。

p_m² > n のとき g(n, m) = f(n, m) であることに注意する。したがって、g を f で初期化し、m を大きい順に小さい順に再帰的に計算すればよい。最終的な答えは g(n, 0) となる。利点は、計算後、g(⌊n/i⌋, 0) のすべての値を含む配列が得られるため、⌊n/i⌋ を求める際に再計算する必要がないことである。

コードの実装は g(n, m-1) = g(n, m) + Σ_c F(p_m^c)(g(⌊n / p_m^c⌋, m) - psum(m)) + F(p_m^c+1) となり、主に列挙の上限を減らすために便利である。

• 複雑度

不明。皆が O(n^0.75 / ln n) と言うが、c を列挙する分、複雑度が増大するかもしれないと感じる。しかし、複雑度が O(n^0.75) を超えないことは証明できる。

応用

最も一般的な応用は、様々な積性関数の和を求める問題である。しかし、基本的にはテンプレートを当てはめるだけである。

min25筛の柔軟な応用を考察する問題もある。

素数カウント

n 以下の素数のうち、mod で r に等しい数がいくつあるかを求める。n ≤ 10¹⁰, r ≤ 13

min25筛のステップ1の考察。

最初のDP配列に次元を一つ追加し、f(n, m, r) を、2 から n までの数のうち、前 m 個の素数を含まないか、素数自体である数の個数と定義する。遷移は p_m を取り除くことで行う：

次大素因数の和

sum_{i ≤ n} f(i) を求める。ここで f(i) は i の非厳密な次大素因数（最大の素因数を除いた後の最大の素因数）を表す。なければ 0 とする。n ≤ 10¹¹。

min25筛のステップ2の考察。

min25筛のステップ2は、本質的に各数の各素因数を大きい順にふるい落としている。最大の素因数は個別に処理される。次大素因数をふるい落とすには、現在考慮している因数が次大素因数であると指定し（そうでなければ継承）、それより大きい数は単一の素数しかないため、高速な素数カウントが可能となる。

コード

min25筛

// divcntk
int sqrt_n, id_count, id_map[N];
ull id[N], n, power_k, dp[N];

inline int get_index(ull x) {
    return id_map[x <= sqrt_n ? x : sqrt_n + n / x];
}

inline void precompute_dp() {
    for (ull l = 1, r; l <= n; l = r + 1) {
        ull zhi = n / l;
        r = n / zhi;
        id[++id_count] = zhi;
        id_map[zhi <= sqrt_n ? zhi : sqrt_n + n / zhi] = id_count;
        dp[id_count] = zhi - 1;
    }
    for (int j = 1; 1ull * prime[j] * prime[j] <= n; ++j) {
        for (int i = 1; id[i] >= 1ull * prime[j] * prime[j]; ++i) {
            dp[i] = dp[i] - dp[get_index(id[i] / prime[j])] + j - 1;
        }
    }
}

ull calculate_g(ull n, int j) {
    if (prime[j] >= n) return 0;
    ull result = (dp[get_index(n)] - j) * (power_k + 1);
    for (int k = j + 1; 1ull * prime[k] * prime[k] <= n; ++k) {
        ull tmp = prime[k];
        for (int e = 1; tmp <= n; ++e, tmp *= prime[k]) {
            result += (power_k * e + 1) * (calculate_g(n / tmp, k) + (e > 1));
        }
    }
    return result;
}

inline void solve() {
    read(n), read(power_k);
    sqrt_n = sqrt(n);
    precompute_dp();
    ull result = calculate_g(n, 0) + 1;
    printf("%llu
", result);
}

洲閣筛

// divcnt3
inline void sieve_primes(const int up) // 素数をふるいにかける
inline int get_index(ll x) { return x <= sqrt_n ? x : tot - n / x + 1; }
ll range[N], tot, dp[N];

inline ll calculate_value(ll k) { return 3 * k + 1; } // h(p^k)

inline void solve() {
    read(n);
    sqrt_n = sqrt(n);
    sieve_primes(sqrt_n + 200);
    for (ll l = 1; l <= n; l = range[tot] + 1)
        range[++tot] = n / (n / l), dp[tot] = range[tot] - 1;
    int ptr = 1;
    for (int j = 1; prime[j] * prime[j] <= n; ++j) {
        while (prime[j] * prime[j] > range[ptr]) ++ptr;
        for (int i = tot; i >= ptr; --i) {
            dp[i] -= dp[get_index(range[i] / prime[j])] - j + 1;
        }
    }
    for (int i = 1; i <= tot; ++i) dp[i] *= 4;
    ptr = tot + 1;
    for (int j = prime_count; j; --j) {
        while (prime[j] * prime[j] <= range[ptr - 1]) --ptr;
        for (int i = tot; i >= ptr; --i) {
            ll nw = prime[j], nww = nw * nw;
            for (int c = 1; nww <= range[i]; ++c, nw = nww, nww *= prime[j]) {
                dp[i] += calculate_value(c) * (dp[get_index(range[i] / nw)] - 4 * j) + calculate_value(c + 1);
            }
        }
    }
    ll result = dp[tot] + 1;
    printf("%lld
", result);
}

注意：

i と r[i] を区別すること！