The 2021 CCPC Guangzhou Onsite - K. Magus Night 题解

题意

Marisa Kirisame 有 \(n\) 个魔法珠，他可以进行一次充能，使得每个珠子 \(s_i\) 随机充能到 \([1,m]\) 的一个值，如果这些珠子的魔法值的 \(gcd\leq p\) 且 \(lcm\geq q\) 那么这算一次成功的充能，带来 \(\prod_{i=0}^n s_i\) 的贡献，否则贡献为 \(0\) ，记贡献期望值为 \(E\) ，要求输出 \(E \cdot m^n \ mod\ \ 998\ 244\ 353\) 的结果

思路

容斥套容斥套容斥，三层容斥，赛上没想到这么多层，也没想到gcd限制+lcm限制下是可以在2s内跑出来的，麻了

---

划分

先来个容斥：由于 \(lcm\geq p\) 的条件不好求贡献，考虑转换为 \(lcm<p\) 后枚举 \(lcm\)（因为枚举 \(lcm\geq p\) 的话，要注意到 \(lcm\) 会超过 \(m\) ）

记（写法可能不标准）

对数组 \(S = \{s_1,s_2,s_3, \dots ,s_n\}\)

\[\begin{aligned} Z &= \{S\ |\ 1\leq s_i\leq m\} \\ A &= \{Z\ |\ gcd(S) \leq q\ \land lcm(S)\geq p\} \\ B &= \{Z\ |\ gcd(S) > q \} \\ C &= \{Z\ |\ gcd(S) \leq q\ \land lcm(S)<p\} \end{aligned} \]

故

\[A = Z - B - C \]

Z

\[\begin{aligned} &(1+2+3+\dots +m)^n \\ =&\big(\frac{1}{2}\cdot m \cdot (m+1)\big)^n \end{aligned} \]

复杂度 \(O(log\ n)\)

B

如果你不会求这个……那么你应该学习 luogu P4450 后再来看这篇题解

总体来说，就是枚举 \(s_i\) 为 \(1,2,3,\dots ,m\) 的倍数后乘上一个对应的莫比乌斯函数

贡献即为（ \(k\) 表示 \(gcd=k\) 时的贡献）

\[\begin{aligned} & \sum_{k=q+1}^m \sum _{i=1} ^{\lfloor \frac{m}{k} \rfloor} \mu(i) \cdot (i+2i+3i+ \dots + {\lfloor \frac{m}{i\cdot k} \rfloor}\cdot i)^n \\ =& \sum_{k=q+1}^m \sum _{i=1} ^{\lfloor \frac{m}{k} \rfloor} \mu(i) \cdot \big(\frac{1}{2}\cdot ik \cdot {\lfloor \frac{m}{i\cdot k} \rfloor} ({\lfloor \frac{m}{i\cdot k} \rfloor}+1)\big) ^ n \end{aligned} \]

复杂度\(O(m\cdot log\ m \cdot log\ n)\)

C

前方容斥套容斥警告

记 \(f(j)\) 表示 \(lcm(S)=j\) 时的贡献

\[\begin{aligned} f(j) &= \sum_{d|j}\mu(\frac{j}{d}) \cdot \Big(\sum d\Big)^n \\ \end{aligned} \]

后面那个 \(\sum d\) 可以预处理一下

所以 \(C\) 的贡献是

\[\begin{aligned} & \sum _{k=1}^{m} \sum _{i=1}^{\lfloor\frac{m}{k} \rfloor} \mu(i) \cdot \sum _{j=1}^{\lfloor \frac{p-1}{i\cdot k}\rfloor} f(j) \cdot (i\cdot k)^n \\ =& \sum _{k=1}^{m} \sum _{i=1}^{\lfloor\frac{m}{k} \rfloor} \mu(i) \cdot (i\cdot k)^n \cdot \sum _{j=1}^{\lfloor \frac{p-1}{i\cdot k}\rfloor} \sum_{d|j}\mu(\frac{j}{d}) \cdot \Big(\sum d\Big)^n \end{aligned} \]

复杂度O(能过)

复杂度我也分析不清楚，但是不会超过\(O(m\cdot log\ m \cdot log^2\ p)\) ，预处理 \(\sum d\) 下

代码

#include<bits/stdc++.h>
using namespace std;
using i64 = long long;

const int N = 2e5+6;
vector<int> divs[N];
i64 dsum[N] = {0};
int prms[N],cnt=0,mob[N];
bool notP[N] = {false};
const int MOD = 998244353;

i64 qpow(i64 b,int p) {
    i64 ret=1;
    for(;p;p>>=1,b=b*b%MOD) if(p&1) ret = ret*b%MOD;
    return ret;
}

void get_mob() {
    mob[1] = 1;
    for(int i=2;i<N;++i) {
        if(!notP[i]) {
            prms[cnt] = i;
            ++cnt;
            mob[i] = -1;
        }

        for(int j=0;j<cnt;++j) {
            if(i*prms[j]>=N) break;

            mob[i*prms[j]] = -mob[i];
            notP[i*prms[j]] = true;

            if(i%prms[j]==0) {
                mob[i*prms[j]] = 0;
                break;
            }
        }
    }
}

void norm(i64 &x) {
    while(x >= MOD) x -= MOD;
    while(x < 0) x += MOD;
}

void get_d() {
    for(int i=1;i<N;++i) {
        for(int j=i;j<N;j+=i) {
            divs[j].push_back(i);
            dsum[j] += i;
            norm(dsum[j]);
        }
    }
}

int main(int argc, char const *argv[])
{
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr); cout.tie(nullptr);

    get_mob();
    get_d();

    i64 n;
    int m,p,q;
    cin >> n >> m >> p >> q;

    for(int i=0;i<N;++i) {
        dsum[i] = qpow(dsum[i],n);
    }

    i64 ans2=0,ans3=0,ans4=0;

    // 2
    for(int k=q+1;k<=m;++k) {
        int M = m/k;
        for(int i=1;i<=M;++i) {
            ans2 += mob[i] * qpow(1ll*(M/i)*(M/i+1)/2%MOD*i%MOD*k%MOD, n);
            norm(ans2);
        }
    }

    // cout << ans2 << endl;

    // 3
    for(int k=1;k<=q;++k) {
        int M = m/k;
        for(int i=1;i<=M;++i) {
            i64 tmp = 0;
            int P = (p-1)/i/k;
            for(int j=1;j<=P;++j) {
                for(int d : divs[j]) {
                    tmp += mob[j/d] * dsum[d] % MOD;
                    norm(tmp);
                }
            }
            ans3 += mob[i] * tmp * qpow(i*k%MOD,n) % MOD;
            norm(ans3);
        }
    }

    // cout << ans3 << endl;

    // 4
    ans4 = qpow(1ll*m*(m+1)/2%MOD,n);

    i64 ans = ans4 - ans2 - ans3;
    norm(ans);
    cout << ans;

    return 0;
}