CCPC广州2021 K.Magus Night

这题神啊…虽然感觉只是步骤很多而已…是因为这题不会没有金也认了啊Orzzz…

• 题目大意：有长度为n的序列，每个数字在$[1,m]$之间，在$m^n$种序列中，定义"好"的序列为：序列的$lcm>=p,gcd<=q$。一个序列的权值定义为：所有元素的乘积。对于所有"好"的序列，求权值和。
• $n<=998244351,m,p,q<=2e5$

题目要求求$lcm>=p,gcd<=q$的好序列。
所以我们将所有序列分为三类：
$A:gcd>q$
$B:gcd<=q,lcm<p$
$C:gcd<=q,lcm>=p$（目标要求）
由于全集是，所有数组的和，为$(m\ast (m+1)/2{)}^n$，所以只要求出AB集合的权值减去即可。

对于$A$集合来说：
考虑$dpa[i]$表示，gcd恰好为i的序列权值和。
由于$gcd\%i==0$的序列权值和可以快速得到$(i\ast (m/i\ast (m/i+1))/2{)}^n$
所以我们从i开始，枚举i的倍数减去即可(基本容斥)

对于$B$集合来说：
考虑$dpb[i]$表示，$lcm=i,gcd=1$的序列权值和。
那么可以把i分解质因数，对于每一个质因子分别考虑。
可以发现，对于一种质因子$p$，出现了$c$次，一定要存在$p^0,p^c$两种数字才可满足。
所以我们容斥一下，把四种情况的贡献都算出来就行了。
不同质因子之间答案相乘即可。

我们有了$dpb$数组，就可以通过枚举$gcd$，通过枚举倍数来统计B集合，复杂度为调和级数。
最后使用全集-AB集合的权值和即为答案。

代码实现如下：

#pragma GCC optimize("O3")
#pragma GCC target("avx2")
#include<bits/stdc++.h>
#define ll long long
#define maxn 400005
#define inf 1e9
#define ins insert
#define pb push_back
#define vi vector <int>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
#define per(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

inline ll read()
{
	ll x=0,w=1; char c=getchar();
	while(c<'0'||c>'9') {if(c=='-') w=-1; c=getchar();}
	while(c<='9'&&c>='0') {x=(x<<1)+(x<<3)+c-'0'; c=getchar();}
	return w==1?x:-x;
}

ll dpa[maxn],dpb[maxn],P[maxn],n,m,p,q,mn[maxn],ANS;
bool is[maxn];
const ll mod=998244353;

inline ll pw(ll a,ll b)
{
	ll ans=1,base=a;
	while(b)
	{
		if(b&1) ans=(ans*base)%mod;
		base=(base*base)%mod; b>>=1;
	}
	return ans;
}

int main()
{
	n=read(); m=read(); p=read(); q=read();
	rep(i,1,2*m) P[i]=pw(i,n);
	per(i,m,1)
	{
		dpa[i]=pw((m/i*(m/i+1)/2)%mod*i%mod,n);
		for(ll j=2*i;j<=m;j+=i) dpa[i]=(dpa[i]-dpa[j]+mod)%mod;
	}
	rep(i,2,maxn-5)
	{
		if(is[i]) continue; mn[i]=i;
		for(ll j=2*i;j<=maxn-5;j+=i) is[j]=1,mn[j]=i;
	}
	
	rep(i,1,m)
	{
		int tmp=i; ll ans=1;
		while(tmp!=1)
		{
			int x=mn[tmp],ret=0,nw=1;
			while(tmp%x==0)
			{
				ret+=nw;
				tmp/=x;
				nw*=x;
			}
			ll qq=(P[ret+nw]-P[ret+nw-1]-P[ret]+P[ret-1]);
			qq=(qq%mod+mod)%mod; ans=(ans*qq)%mod;
		}
		dpb[i]=ans;
	}
	
	ANS=pw(((m*(m+1))/2)%mod,n);
	rep(i,q+1,m) ANS=(ANS-dpa[i]+mod)%mod;
	rep(i,1,q)
	{
		for(ll j=i,k=1;j<p;j+=i,k++)
			ANS=(ANS-P[i]*dpb[k]%mod+mod)%mod;
	}
	cout<<ANS<<endl;
	return 0;
}