CCPC2021 廣州 K. Magus Night

題意

給定整數區間 \([1,m]\) ，從中可重復的選擇 \(n\) 個數，形成一個數列 \(\{a_n\}\) 。問：所有滿足 \(\gcd(a_1,...,a_n)\le q\) 並且 \(\text{lcm}(a_1,...,a_n)\ge p\) 的數列的乘積和。

題解

官方題解其實已經很明了了，我這里再做個翻譯。題目要求的是 \(\gcd\le q\) 且 \(\text{lcm} \ge p\) 的所有數列的乘積和。根據 \(\gcd|\text{lcm}\) 這一性質可以枚舉 \((\gcd,\text{lcm})\) 數對來計算，但此題的 \(\text{lcm}\) 會很大，無法直接枚舉所有大於 \(p\) 的 \(\text{lcm}\) ，所以需要轉換思路，大的枚舉不完，但小的好枚舉，於是我們做這么一個容斥：

\[(\gcd\le q,\text{lcm}\ge p)=S-(\gcd>q)-(\gcd\le q,\text{lcm}<p) \]

其中 \((a,b)\) 表示滿足 \(a\) 和 \(b\) 條件下的所有數列的乘積和， \(S\) 代表全體數列的乘積和。

全體數列的乘積和可以這么表示：每一個數的取值范圍都是 \([1,m]\) ，利用分配律，可知

\[S=(1+2+...+m)(1+2+...+m)...(1+2+...+m)=(1+2+...+m)^n \]

然后再求 \((\gcd>q)\) 。這是個經典問題，很容易由容斥算出來。

然后就是 \((\gcd\le q,\text{lcm}<p)\) 。考慮上述的分配律，我們也可以把滿足 \(\gcd=g,\text{lcm}=l\) 的所有數列的乘積和表示成一些數乘積的和。考慮到唯一分解定理，我們使用素數來進行表示。對於兩個數 \(a,b\) ，把他們寫成唯一分解的形式：

\[a=p_1^{k_1}p_2^{k_2}...p_s^{k_s}\\ b=p_1^{t_1}p_2^{t_2}...p_s^{t_s}\\ \gcd(a,b)=p_1^{\min(k_1,t_1)}p_2^{\min(k_2,t_2)}...p_s^{\min(k_s,t_s)}\\ \text{lcm}(a,b)=p_1^{\max(k_1,t_1)}p_2^{\max(k_2,t_2)}...p_s^{\max(k_s,t_s)} \]

由此可知，一對 \(\gcd,\text{lcm}\) 便確定了素數冪次的上界和下界，利用上述的分配律計算貢獻即可。計算的時候需要注意的是，上界和下界必須取到，可以考慮下面的容斥來進行計算：

設下界為 \(L\) ，上界為 \(R\) ，那么

\[ans(L,R)=ans(L\le x\le R)-ans(L\le x\le R-1)-ans(L+1\le x\le R)+ans(L+1\le x\le R-1) \]

AC代碼

這份代碼沒有優化，1887ms，過的非常危險。

#include <bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<b;i++)
using namespace std;

using ll = long long;

constexpr int N = 2e5+5, P=998244353;
ll f[N], mind[N];

ll qpow(ll a, ll b){
	ll res=1;
	for(;b;b>>=1,a=a*a%P)if(b&1)res=res*a%P;
	return res;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("1.in", "r", stdin);
#endif // !ONLINE_JUDGE
    cin.tie(nullptr)->ios::sync_with_stdio(false);
	rep(i,2,N)if(!mind[i])for(int j=i;j<N;j+=i)mind[j]=i;
    ll n,m,p,q; cin>>n>>m>>p>>q;
	ll inv2=qpow(2,P-2);
	ll ans=qpow((m*(m+1)%P*inv2)%P,n);
	for(int i=m;i>q;i--){
		ll cnt=m/i;
		f[i]=qpow(i,n)*qpow((cnt*(cnt+1)%P*inv2%P),n)%P;
		for(int j=i+i;j<=m;j+=i)f[i]=(f[i]-f[j]+P)%P;
		ans=(ans-f[i]+P)%P;
	}//cout<<ans<<endl;
	rep(l,1,m+1){
		f[l]=1;
		for(int tmp=l,pr;tmp!=1;){
			pr=mind[tmp];
			ll x=0,y=1;
			for(;tmp%pr==0;x=(x+y)%P,y=y*pr%P)tmp/=pr;
			ll temp=(qpow(x+y,n)-qpow(x,n)-qpow(x+y-1,n)+qpow(x-1,n))%P;
			f[l]=(f[l]*temp)%P;
		}
	}
	rep(i,1,q+1)for(int j=i;j<p;j+=i)ans=(ans-qpow(i,n)*f[j/i]%P+P)%P;
	cout<<ans;
    return 0;
}