杜教篩

数論関数 \(f(n)\) の累積和 \(F(n)=\sum_{i=1}^n f(i)\) を求めるとします。ここで、以下の条件を満たす積性関数 \(g(n)\) が存在する場合:

• \(g(n)\) の累積和 \(G(n)=\sum_{i=1}^n g(i)\) が効率的に計算可能
• ディリクレ積 \(h = f \ast g\) の累積和 \(H(n)=\sum_{i=1}^n h(i)\) が効率的に計算可能

以下の関係式が導出されます:

\[\begin{aligned} \sum_{i=1}^n (f \ast g) &= \sum_{i=1}^n \sum_{d|i} f(d)g\left(\frac{i}{d}\right) \\ &= \sum_{d=1}^n g(d) \sum_{k=1}^{\lfloor n/d \rfloor} f(k) \\ &= F(n) + \sum_{d=2}^n g(d) F\left(\left\lfloor \frac{n}{d} \right\rfloor\right) \end{aligned}\]

したがって:

\[F(n) = H(n) - \sum_{d=2}^n G(d) F\left(\left\lfloor \frac{n}{d} \right\rfloor\right)\]

実装例 (\(f(n) = n \varphi(n)\) の場合):

constexpr int mod = 1000000007;
constexpr int inv6 = 166666668; // 1/6 mod 1000000007

long long calcG(long long n) {
    if (n <= N) return sg[n];
    if (G.find(n) != G.end()) return G[n];
    
    long long res = n % mod;
    res = (res * (n + 1) % mod) % mod;
    res = (res * (2*n + 1) % mod) % mod;
    res = (res * inv6) % mod;
    
    for (long long l = 2, r; l <= n; l = r + 1) {
        long long q = n / l;
        r = n / q;
        long long term = (r - l + 1) % mod;
        term = (term * (l + r) % mod) % mod;
        term = (term * calcG(q) % mod) % mod;
        term = (term * inv2) % mod;
        res = (res - term + mod) % mod;
    }
    return G[n] = res;
}

Powerful Number 篩

積性関数 \(f(n)\) の累積和 \(F(n)\) を求めるとします。以下の条件を満たす積性関数 \(g(n)\) が存在する場合:

• \(g(n)\) の累積和 \(G(n)\) が効率的に計算可能
• 任意の素数 \(p\) に対して \(g(p) = f(p)\)

積性関数 \(h\) を \(f = h \ast g\) と定義すると:

\[\begin{aligned} F(n) &= \sum_{i=1}^n (h \ast g) \\ &= \sum_{d=1}^n h(d) G\left(\left\lfloor \frac{n}{d} \right\rfloor\right) \end{aligned}\]

重要な性質:

• 任意の素数 \(p\) で \(h(p) = 0\)
• Powerful Number の総数は \(O(\sqrt{n})\)

\(h(p^k)\) の計算:

\[f(p^k) = \sum_{i=0}^k g(p^i) h(p^{k-i})\]

再帰的実装例 (\(f(p^k) = p^k(p^k - 1)\)):

void computePN(int idx, long long x, long long h_val) {
    if (idx > prime_cnt || x * primes[idx] > n || x * primes[idx] * primes[idx] > n) return;
    
    computePN(idx + 1, x, h_val);
    
    long long base = primes[idx] * primes[idx];
    x *= base;
    vector<long long> H = {1, 0};
    long long p_power = base % mod;
    
    for (int exp = 2; x <= n; exp++, x *= primes[idx], p_power = (p_power * primes[idx]) % mod) {
        long long val = p_power * (p_power - 1 + mod) % mod;
        long long g_term = 1;
        
        for (int t = 1; t <= exp; t++, g_term = (g_term * primes[idx]) % mod) {
            long long coef = g_term * g_term % mod;
            coef = coef * primes[idx] % mod;
            coef = coef * (primes[idx] - 1) % mod;
            val = (val - coef * H[exp - t] % mod + mod) % mod;
        }
        H.push_back(val);
        result = (result + h_val * val % mod * calcG(n / x) % mod) % mod;
        computePN(idx + 1, x, h_val * val % mod);
    }
}

Min_25 篩

積性関数 \(f(n)\) の累積和 \(F(n)\) を効率的に計算します。状態を以下で定義:

• \(S(n, j)\): \(n\) 以下で最小素因数が \(p_j\) より大きい整数に対する \(f(x)\) の総和

最終的に \(F(n) = S(n, 0)\) となります。計算は素数部分と合成数部分に分割:

1. 素数部分: \(g(n, \infty)\) で表現
2. 合成数部分: 最小素因数 \(p_k\) の冪乗 \(p_k^e\) について計算

実装例:

long long computeS(long long n, int j) {
    if (n < primes[j]) return 0;
    int idx = index_map[n];
    long long prime_sum = (g2[idx] - g1[idx] - (sp2[j] - sp1[j]) + 2*mod) % mod;
    
    for (int k = j + 1; k <= prime_cnt && (long long)primes[k]*primes[k] <= n; k++) {
        for (long long e = primes[k]; e <= n; e *= primes[k]) {
            long long term = e % mod;
            term = term * (e % mod - 1 + mod) % mod;
            prime_sum = (prime_sum + term * (computeS(n/e, k) + (e != primes[k]))) % mod;
        }
    }
    return prime_sum;
}

int main() {
    for (long long l = 1, r; l <= n; l = r + 1) {
        long long x = n / l;
        r = n / x;
        values[++idx_cnt] = x;
        index_map[x] = idx_cnt;
        
        long long term1 = (x % mod) * ((x + 1) % mod) % mod;
        term1 = term1 * inv2 % mod;
        g1[idx_cnt] = (term1 - 1 + mod) % mod;
        
        long long term2 = (x % mod) * ((x + 1) % mod) % mod;
        term2 = term2 * ((2*x + 1) % mod) % mod;
        term2 = term2 * inv6 % mod;
        g2[idx_cnt] = (term2 - 1 + mod) % mod;
    }
    
    for (int i = 1; i <= prime_cnt; i++) {
        long long p_sq = (long long)primes[i] * primes[i];
        for (int j = 1; j <= idx_cnt && p_sq <= values[j]; j++) {
            long long q_val = values[j] / primes[i];
            int q_idx = index_map[q_val];
            
            g1[j] = (g1[j] - primes[i] * (g1[q_idx] - sp1[i-1] + mod) % mod + mod) % mod;
            g2[j] = (g2[j] - (long long)primes[i]*primes[i] % mod * (g2[q_idx] - sp2[i-1] + mod) % mod + mod) % mod;
        }
    }
    cout << (computeS(n, 0) + 1) % mod << endl;
}