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

title: The 2021 CCPC Guangzhou Onsite - K. Magus Night 題解
date: 2021-11-18 18:19:00
tags: [inclusion-exclusion, number theory, math, team training]
mathjax: true

題意

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;
}