5663. 【GDOI2018Day1模拟4.17】呼吸决定（杜教筛 + 拉格朗日插值法）

佚名 6年前 (2019-04-12) 随笔 896人围观抢沙发百度已收录

Problem

\[\sum_{i=1}^n\mu(i)i^m\]

\(n\le 10^9, m\le 10^5\)

SRE实战互联网时代守护先锋，助力企业售后服务体系运筹帷幄！一键直达领取阿里云限量特价优惠。

Solution

直接上杜教筛就好了。

如果让\(f(i)=\mu(i)i^m,g(i)=i^m\)，可以得到\(f·g=h=[n=1]\)。

然后常规套路：\[S(n)=\sum_{i=1}^nf(i)\]，然后会发现变成一个\[S(n)=1-\sum_{d=2}^nd^mS(\frac{n}{d})\]

套一个拉格朗日插值就好了。

关键是卡空间。。

这个需要强大的卡卡能力，然后注意些数组可以省去，然后尽量的卡一卡就好了。

Code

#include <bits/stdc++.h>

#define F(i, a, b) for (int i = (a); i <= (b); i ++)
#define G(i, a, b) for (int i = (a); i >= (b); i --)

const int Mo = 998244353, N = 5e6 + 10, K = 1e5 + 10, T = 1000;

int n, n2, k, m, ny[K], z[N / 13];
int mu[N / 5], y[N], S[N / 5], vis[T], pre[K], suf[K], Sum[T / 5];

using namespace std;

int ksm(int x, int y) {
    int ans = 1;
    for (; y ; y >>= 1, x = (1ll * x * x) % Mo)
        if (y & 1)
            ans = (1ll * ans * x) % Mo;
    return ans;
}

int Doit(int x) {
    if (x == 1) return 1;
    if (x < N / 5) return S[x];
    int num = n / x;
    if (vis[num])
        return vis[num];
    int Ans = 0;
    for (int i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        int xx = n / i, yy = n / (j + 1);
        int Sx = (i - 1 < N) ? y[i - 1] : Sum[xx], Sy = (j < N) ? y[j] : Sum[yy];
        Ans = (Ans + 1ll * (Sy - Sx) * Doit(x / i)) % Mo;
    }
    return vis[num] = 1 - Ans;
}

int Get(int n) {
    int sum = 0;
    pre[0] = suf[m + 1] = 1;
    F(j, 1, m)
        pre[j] = 1ll * pre[j - 1] * (n - j) % Mo;
    G(j, m, 1)
        suf[j] = 1ll * suf[j + 1] * (n - j) % Mo;
    F(i, 1, m)
        if ((m - i) & 1)
            sum = (sum - 1ll * y[i] * pre[i - 1] % Mo * suf[i + 1] % Mo * ny[i - 1] % Mo * ny[m - i] % Mo) % Mo;
        else
            sum = (sum + 1ll * y[i] * pre[i - 1] % Mo * suf[i + 1] % Mo * ny[i - 1] % Mo * ny[m - i] % Mo) % Mo;
    (sum += Mo) %= Mo;
    return sum;
}

int main() {
    freopen("calc.in", "r", stdin);
    freopen("calc.out", "w", stdout);
    
    scanf("%d%d", &n, &k), n2 = int(sqrt(n));
    F(i, 2, N - 1) {
        if (!y[i]) {
            z[++ z[0]] = i, y[i] = ksm(i, k);
            if (i < N / 5) mu[i] = - 1;
        }
        F(j, 1, z[0]) {
            long long x = 1ll * z[j] * i;
            if (x >= N) break;
            y[x] = (1ll * y[z[j]] * y[i]) % Mo;
            if (i % z[j] == 0) break;
            if (x < N / 5)
                mu[x] = mu[z[j]] * mu[i];
        }
    }
    

    y[1] = S[1] = 1;
    F(i, 2, N / 5 - 1)
        S[i] = (long long) (S[i - 1] + 1ll * mu[i] * y[i]) % Mo;
    F(i, 2, N - 1)
        y[i] = (long long) (y[i] + y[i - 1]) % Mo;

    ny[0] = 1, m = k + 2;
    F(i, 1, m)
        ny[i] = (1ll * ny[i - 1] * i) % Mo;
    ny[m] = ksm(ny[m], Mo - 2); 
    G(i, m - 1, 1)
        ny[i] = (1ll * ny[i + 1] * (i + 1)) % Mo;

    for (int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        if (i >= N)
            Sum[n / i] = Get(i) - ksm(i, k);
    }
    Sum[0] = Get(n + 1) - ksm(n + 1, k);

    printf("%d\n", (Doit(n) + Mo) %Mo);
}