CF1559E -Mocha and Stars (容斥+dp+莫比乌斯函数)

题目

给定n个区间\([l_i,r_i]\)和数\(m\)，问有多少组\((a_1,a_2,...,a_i)\)满足：

• \(a_i\in [l_i,r_i]\)

• \(\sum\limits_{i=1}^na_i\le m\)

• \(\gcd(a_1,...,a_n)=1\)

数据范围：\(n\le 50\)，\(m \le 10^5\)

题解

先想想不考虑\(\gcd\)的限制的情况下怎么做。就是一个简单的$$dp$$，$$dp[i][j]$$代表前\(i\)个区间总和为\(j\)的方案数。

\[dp[i][j]=\sum\limits_{k=j-r_i}^{j-l_i}{dp[i-1][k]} \]

时间复杂度为\(O(n\cdot m)\)，\(m\)为最大值。

考虑枚举\(k|\gcd(a_1,...,a_n)\)的\(k\)的方案数，可以容斥，容斥因子是\(\mu\)函数。对于\(k\)，最大值变为\(\frac{m}{k}\)，总时间复杂度为\(O(n\cdot m \cdot \log m)\)。

这个容斥因子是\(\mu\)函数的似乎是利用本质不同的质因数个数按照奇减偶加原则容斥。

#include <bits/stdc++.h>

#define endl '\n'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const int N = 3e5 + 10;
const int M = 998244353;
const double eps = 1e-5;
int rev[N];

typedef pair<int, int> PII;
PII seg[N];
int cnt;
int prime[N];
int isnp[N];
int mu[N];
ll dp[60][N];

ll solve(int p, int n, int m) {
    int mx = m / p;
    for(int i = 0; i <= mx; i++) {
        dp[0][i] = 1;
    }
    for(int i = 1; i <= n; i++) {
        bool flag = true;
        int l = seg[i].first / p + (seg[i].first % p != 0), r = seg[i].second / p;
        if(l > r) return 0;
        for(int j = 1; j <= mx; j++) {
            dp[i][j] = ((j - l < 0 ? 0 : dp[i - 1][j - l]) - (j - r - 1 < 0 ? 0 : dp[i - 1][j - r - 1]) + M) % M;
            dp[i][j] = (dp[i][j] + dp[i][j - 1]) % M;
        }
    }
    return dp[n][mx];
}


int main() {
    int n, m;
    mu[1] = 1;
    isnp[1] = 1;
    prime[0] = 1;
    for(int i = 2; i < N; i++) {
        if(!isnp[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; j <= cnt && i * prime[j] < N; j++) {
            isnp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    cin >> n >> m;
    for(int i = 1; i <= n; i++) {
        cin >> seg[i].first >> seg[i].second;
    }
    ll ans = 0;
    for(int i = 1; i <= m; i++) {
        if(!mu[i]) continue;
        ans += mu[i] * solve(i, n, m);
        ans %= M;
    }
    cout << (ans + M) % M << endl;
}