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}\) 便确定了素数幂次的上界和下界，利用上述的分配律计算贡献即可。

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