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\)为最大值。

然后就是经典的莫比乌斯反演

\[\begin{aligned} ans&=\sum\limits_{x_1=1}^{[l_1,r_1]}{\sum\limits_{x_2=1}^{[l_2,r_2]}{...\sum\limits_{x_n=1}^{[l_n,r_n]}{[x_1+...+x_n\leq m]\cdot [\gcd(x_1,...,x_n)]}}} \\ &=\sum\limits_{x_1=1}^{[l_1,r_1]}{\sum\limits_{x_2=1}^{[l_2,r_2]}{...\sum\limits_{x_n=1}^{[l_n,r_n]}{[x_1+...+x_n\leq m]\cdot \sum\limits_{d|x_1,d|x_2,...,d|x_n}{\mu(d)}}}} \\ &=\sum_{d=1}^{m}{\mu(d)(\sum\limits_{x_1=1}^{[\frac{l_1}{d},\frac{r_1}{d}]}{\sum\limits_{x_2=1}^{[\frac{l_2}{d},\frac{r_2}{d}]}{...\sum\limits_{x_n=1}^{[\frac{l_n}{d},\frac{r_n}{d}]}{[x_1+...+x_n\leq \frac{m}{d}]}}})} \end{aligned} \]

后面括号那一坨就是前面的\(dp\)可以解决。

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