LOJ6519 魔力环

题面：LOJ

解析

因为是等价环计数，考虑Burnside引理。

设f(i)表示将环分做\(n/i\)个循环的不动点个数。
发现对于\(f(i)\)，其循环长度为\(i\)，那么一定有\(i|m\)，
即：\(i|gcd(n,m)\)，否则没有贡献，所以：

\[Ans=\frac{\sum_{d|gcd(n,m)} f(d) \varphi(\frac{n}{d})}{n} \]

那么如何计算\(f(i)\)呢？此处为方便，特判\(n=m\)的情况，
然后问题就等价于\(i\)个小球成环，给其中\(\frac{mi}{n}\)个小球染色，连续不超过\(k\)个的方案数。

形式化地，令小球个数为\(x\)，黑球个数为\(y\)，答案为\(F(x,y)\)。
注意到这里的环就没有等价问题了，所以考虑破环成链。

现在就是在\(x-y\)个白球里插入\(y\)个黑球，每个空隙不超过\(k\)个，
其中第一个白球前面与最后一个白球后面黑球之和亦不能超过\(k\)个。
现在枚举第一个白球与最后一个白球之间的黑球数量\(i\)，有:

\[F(x,y)=\sum_{i=0}^{k} (i+1)T(y-i,x-y-1) \]

其中，因为环是有序号的，所以要枚举\(i\)个黑球的分配方式，有(i+1)种。

而\(T(x,y)\)则表示将\(x\)个球，分配在\(y\)个盒子里（可以为空），每个盒子不超过\(k\)个的方案数。
现在考虑计算\(T(x,y)\)，用容斥，每次强制\(i\)个盒子超出容量，那么有：

\[T(x,y)=\sum_{i=0}^{min(x/(k+1),y)} (-1)^i {y \choose i} H(x-i*(k+1),y) \]

其中，\(H(x,y)\)表示将\(x\)个相同的小球放入\(y\)个不同的盒子里（可以为空）的方案数，有：

\[H(x,y)={x+y-1 \choose y-1} \]

计算\(T(x,y)\)的复杂度是\(O(x/k)\)，计算\(F(x,y)\)的复杂度是\(O(y)\)，还要枚举约数，
所以计算\(Ans\)的复杂度就是\(O(n+\sigma(n))\)，至此，问题圆满解决。

代码

#include<cstdio>
#define N 100005

using namespace std;

inline int In(){
	char c=getchar(); int x=0;
	for(;c<'0'||c>'9';c=getchar());
	for(;c>='0'&&c<='9';c=getchar()) x=x*10+c-'0';
	return x;
}

const int P=998244353;

inline int min(int a,int b){
	return a<b?a:b;
}

inline int gcd(int a,int b){
	return (b==0)?a:gcd(b,a%b);
}

inline int power(int x,int k){
	if(!x) return 0;
	int s=1,t=x;
	for(;k;k>>=1,t=1ll*t*t%P)
	if(k&1) s=1ll*s*t%P;
	return s;
}

int n,m,K,Gcd,ans,phi[N],fac[N],inv[N],p[N],p_cnt=0; bool vis[N];

inline void Get_phi(){
	phi[1]=1;
    for(int i=2;i<=n;++i){
        if(!vis[i]) p[++p_cnt]=i,phi[i]=i-1;
        for(int j=1;j<=p_cnt;++j){
            if(i*p[j]>n) break; vis[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
}

inline void Get_fac(){
	fac[0]=fac[1]=1;
	for(int i=2;i<=n;++i) fac[i]=1ll*fac[i-1]*i%P;
	inv[n]=power(fac[n],P-2);
	for(int i=n-1;~i;--i) inv[i]=1ll*inv[i+1]*(i+1)%P;
}

inline int C(int x,int y){ return 1ll*fac[x]*inv[y]%P*inv[x-y]%P; }
inline int H(int x,int y){ return C(x+y-1,y-1); }
inline int T(int x,int y){
	int res=0,lim=min(x/(K+1),y);
	for(int i=0,t=1;i<=lim;++i,t=P-t)
	res=(res+1ll*t*C(y,i)%P*H(x-i*(K+1),y)%P)%P;
	return res;
}
inline int F(int x,int y){
	if(x-y==1) return K>=(y-1)?x:0;
	int res=0,lim=min(y,K);
	for(int i=0;i<=lim;++i)
	res=(res+1ll*(i+1)*T(y-i,x-y-1)%P)%P;
	return res;
}

int main(){
	n=In(); m=In(); K=In(); ans=0;
	if(K==1&&m>n/2){ printf("%d\n",0); return 0;}
	if(n==m){ printf("%d\n",(K>=n)); return 0; } 
	Gcd=gcd(n,m); Get_phi(); Get_fac();
	for(int i=1;i<=Gcd;++i) if(Gcd%i==0)
	ans=(ans+1ll*F(n/i,m/i)*phi[i]%P)%P;
	printf("%lld\n",1ll*ans*power(n,P-2)%P);
	return 0;
}