牛客挑战赛46 E 反演

反演

https://ac.nowcoder.com/acm/contest/9510/E

考虑使用min25。
令f(p^c)是当p为质数时的答案,根据积性函数的性质,可以进行min25筛。
但是选择的质数小于等于m时,f(p^c) != c + 1。
处理这些奇异值即可进行min25筛。

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

typedef long long ll;

const int N = 1e6 + 100;
const int MOD = 998244353;

int Sqr, K, m, tot, id1[N], id2[N];
bool is_prime[N];
ll n, pri[N], w[N], g[N], spf[N];
ll tt[110], cc[110], rev[110], R;

void Sieve(int n) {
    memset(is_prime, true, sizeof(is_prime));
    is_prime[1] = false;
    for (int i = 2; i <= n; i++) {
        if (is_prime[i]) {
            pri[++tot] = i;
            spf[tot] = spf[tot - 1] + 1;
        }
        for (int j = 1; i * pri[j] <= n; j++) {
            is_prime[i * pri[j]] = false;
            if (i % pri[j] == 0) break;
        }
    }
}

ll qpow(ll x, ll n) {
    ll res = 1;
    while (n > 0) {
        if (n & 1) res = res * x % MOD;
        x = x * x % MOD;
        n /= 2;
    }
    return res;
}

ll F(ll p, ll e) {
    if (p <= K) return rev[p] * ((e + 1) * (tt[p] + 1) % MOD + cc[p]) % MOD; //奇异值
    return e + 1;
}

ll S(ll x, int y) {
    if (x <= 1 || pri[y] > x) return 0;
    int k = (x <= Sqr) ? id1[x] : id2[n / x];
    ll res = g[k] - spf[y - 1];
    for (int i = y; i <= tot && 1LL * pri[i] * pri[i] <= x; i++) {
        ll p1 = pri[i], p2 = pri[i] * pri[i];
        for (int e = 1; p2 <= x; e++, p1 = p2, p2 *= pri[i]) {
            res = (res + S(x / p1, i + 1) * F(pri[i], e) + R * F(pri[i], e + 1)) % MOD;
        }
    }
    return res;
}

int solve() {
    m = 0; tot = 0;
    Sqr = sqrt(n); Sieve(Sqr);
    for (ll i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i); w[++m] = n / i;
        if (w[m] <= Sqr) id1[w[m]] = m;
        else id2[n / w[m]] = m;
        g[m] = (w[m] - 1) % MOD;
    }
    for (int j = 1; j <= tot; j++) {
        for (int i = 1; i <= m && 1LL * pri[j] * pri[j] <= w[i]; i++) {
            int k = (w[i] / pri[j] <= Sqr) ? id1[w[i] / pri[j]] : id2[n / (w[i] / pri[j])];
            g[i] = (g[i] - (g[k] - spf[j - 1])) % MOD;
        }
    }
        //上半部分是min25模板,下半部分处理奇异值。
    R = 1;
    for (int j = 2; j <= K; j++) {
        if (!is_prime[j]) continue;
        for (int i = j * 2; i <= K; i += j) is_prime[i] = false;
        for (int i = 2; i <= K; i++) {
            int x = i;
            while (x % j == 0) {
                x /= j;
                tt[j]++;
                cc[j] += tt[j];
            }
        }
        R = R * (tt[j] + 1 + cc[j]) % MOD;
        rev[j] = qpow(tt[j] + 1 + cc[j], MOD - 2);
    }
    ll sum = 0, cnt = 0;
    for (int i = m, j = 1; i >= 1; i--) {
        while (j + 1 <= K && j + 1 <= w[i]) {
            j++;
            if (!is_prime[j]) continue;
            sum = (sum + R * rev[j] % MOD
                * (2 * (tt[j] + 1) + cc[j])) % MOD;
            cnt++;
        }
        g[i] = (2 * (g[i] - cnt) * R + sum) % MOD;
    }
    for (int j = 1; j <= tot; j++) {
        if (pri[j] <= K) spf[j] = spf[j - 1] + R * rev[pri[j]] % MOD
            * (2 * (tt[pri[j]] + 1) + cc[pri[j]]) % MOD;
        else spf[j] = (spf[j - 1] + R * 2) % MOD;
    }
    return S(n, 1) + R;
}

int main() {
    //freopen("0.txt", "r", stdin);
    scanf("%lld%d", &n, &K);
    ll ans = solve();
    ans = (ans % MOD + MOD) % MOD;
    printf("%lld\n", ans);
    return 0;
}
全部评论

相关推荐

冲芭芭拉鸭:你这图还挺新,偷了。
投递美团等公司10个岗位
点赞 评论 收藏
分享
10-09 22:05
666 C++
找到工作就狠狠玩CSGO:报联合国演讲,报电子烟设计与制造
点赞 评论 收藏
分享
3 收藏 评论
分享
牛客网
牛客企业服务