[CQOI2021模拟]此罪当珠

题目

题目背景
猪八戒调戏嫦娥,玉皇大帝震怒,说:“此罪当诛!”

——然后他就投胎当了猪。

题目描述
有一个 n n n 个珠子的手环,其中 m m m 个是金色的,其余是白色的。

为了产生交错的美感,要求金色的连续段长度不超过 k k k 。即,任意 k + 1 k+1 k+1 个连续的珠子,都不全是金色。

现在请问,不同的手环有多少种?补充一句,手环是可以旋转的,但是不可以翻转;珠子只有颜色的区别。形式化地说,记金色为 1 1 1,白色为 0 0 0,两个手环形成的颜色序列(从任意一个珠子开始,顺时针写下每个珠子的颜色)分别是 ⟨ a 0 , a 1 , … , a n − 1 ⟩ ,    ⟨ b 0 , b 1 , … , b n − 1 ⟩ \langle a_0,a_1,\dots,a_{n-1}\rangle,\;\langle b_0,b_1,\dots,b_{n-1}\rangle a0,a1,,an1,b0,b1,,bn1,如果存在一个正整数 t t t 使得 a ( i + t )   m o d   n = b i a_{(i+t)\bmod n}=b_i a(i+t)modn=bi,那么这两个手环是相同的。

输出方案数对 998244353 998244353 998244353 取模的结构。

数据范围与提示
0 ≤ k ≤ m ≤ n ≤ 1 0 6 0\le k\le m\le n\le 10^6 0kmn106,最多 5 5 5 组数据。

思路

显然是 b u r n s i d e \tt burnside burnside 嘛。我其实也忘了许多,但是记得证明的思路,现场重推了一下。

顺时针旋转 x x x 的操作下,不动点是 a ( i + x )   m o d   n = a i a_{(i+x)\bmod n}=a_i a(i+x)modn=ai,显然就是周期 T = gcd ⁡ ( x , n ) T=\gcd(x,n) T=gcd(x,n) 嘛。只需要确定前 T T T 个珠子的颜色,剩下的就都确定了。

不妨设 m < n m<n m<n,那么至少有一个白色珠子,所以 T T T 中也含有一个白色珠子。真正的手环是 T T T 复制若干次得到的,不难发现,其金色连续段要么完全在 T T T 中,要么是 T T T 首尾相接。

问题转化为,对于一个 T T T 个点的环,求金色连续段长度不超过 k k k 的方案数。不过这个环是不能旋转的(即,旋转相同视为不同方案)。

如果用 d p \tt dp dp 来做,大概就是 f ( i , j ) = ∑ t = 0 k f ( i − 1 , j − t − 1 ) f(i,j)=\sum_{t=0}^{k}f(i-1,j-t-1) f(i,j)=t=0kf(i1,jt1) 的样子。然后你发现它很像母函数。事实上的确如此——将一个金色连续段与它右边的一个白色珠子看成一个物品,每个物品的选择数量不同,则方案不同。用 x x x 的指数代表 金色珠子的数量,则
( ∑ i = 0 k x i ) n − m − 1 \left(\sum_{i=0}^{k}x^i\right)^{n-m-1} (i=0kxi)nm1

为什么是 n − m − 1 n-m-1 nm1 次方?因为一共有 n − m n-m nm 个白色珠子,而最后一个是首尾相接,特殊处理,所以剩下的是 n − m − 1 n-m-1 nm1 次方。(毕竟每个物品恰好提供一个白色珠子。)

首尾相接其实很好处理。对于长度为 i + 1 i+1 i+1(即,有 i i i 的情况,可以在末尾出现 [ 1 , i + 1 ] [1,i+1] [1,i+1] 个珠子(剩下的就放在开头),所以方案数就是 i + 1 i+1 i+1 了。即
G ( x ) = ∑ i = 0 k ( i + 1 ) x i G(x)=\sum_{i=0}^{k}(i+1)x^i G(x)=i=0k(i+1)xi

二者乘起来的第 m m m 项系数就是答案。

那么这里有一个推式子的做法,把 G ( x ) G(x) G(x) 解出来,然后 ∑ i = 0 k x i = 1 − x k + 1 1 − x \sum_{i=0}^{k}x^i=\frac{1-x^{k+1}}{1-x} i=0kxi=1x1xk+1 n − m − 1 n-m-1 nm1 次幂就分别计算分子分母。化简式子之后,大概是 ( ( ( 三项的多项式 ) ) ) × \times × ( ( ( 二项式展开式 ) ) ) × \times × ( ( ( 二项式展开式 ) ) ) 。枚举三项多项式中的一项,再 m k = ⌊ m k ⌋ m_k=\lfloor\frac{m}{k}\rfloor mk=km 枚举 1 − x k + 1 1-x^{k+1} 1xk+1 的二项式展开式中的一项,直接就找到了 1 − x 1-x 1x 的二项式展开式中的一项。

最神奇的是,暴力计算非常快。因为每一项要么取 1 1 − x \frac{1}{1-x} 1x1,要么取 − x k + 1 1 − x \frac{-x^{k+1}}{1-x} 1xxk+1,直接枚举第二类的数量(其不超过 m k m_k mk 个),除以了 x k + 1 x^{k+1} xk+1 后全部变为 1 1 − x {1\over 1-x} 1x1,可以组合数直接计算。——从现实意义来看,就是 容斥原理 计算不超过 k k k 个的方案数。由于 G ( x ) G(x) G(x) 只有 k k k 次,我们只需要求出 ( m − k ) (m-k) (mk) 次到 m m m 次项,这是 O ( k ⋅ m k ) = O ( m ) \mathcal O(k\cdot m_k)=\mathcal O(m) O(kmk)=O(m) 的。

然后枚举 G ( x ) G(x) G(x) 中的一项,乘上上面算出的一项。我们就已经得到了答案。

但是为何我不敢这么打?因为我不会算复杂度。每次计算是 O ( m ) \mathcal O(m) O(m) 的,但是这里的 m m m 不是题目中直接输入的 m m m,而是 T T T 里面有多少个金色的珠子。如果枚举一个 d = n / T d=n/T d=n/T(即周期的个数),可知复杂度为
∑ d ∣ n ,    d ∣ m m d \sum_{d\mid n,\;d\mid m}\frac{m}{d} dn,dmdm

糟糕透顶的情况是 m ∣ n m\mid n mn,此时复杂度是 m m m 的因数和 σ 1 ( m ) = ∏ p i t i + 1 − 1 p i − 1 \sigma_1(m)=\prod{p_i^{t_i+1}-1\over p_i-1} σ1(m)=pi1piti+11,其中 m = ∏ p i t i m=\prod p_i^{t_i} m=piti 表示质因数分解。这东西是近乎线性的(打表证明)。

代码

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
using namespace std;
typedef long long int_;
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int Mod = 998244353;
inline int qkpow(int_ b,int_ q){
	int_ a = 1;
	for(; q; q>>=1,b=b*b%Mod)
		if(q&1) a = a*b%Mod;
	return a;
}

int n, m, k;

namespace MyMath{
	const int MaxN = 3000005;
	int jc[MaxN], inv[MaxN];
	void importC(){
		jc[1] = inv[1] = 1;
		for(int i=2; i<MaxN; ++i){
			jc[i] = 1ll*jc[i-1]*i%Mod;
			inv[i] = inv[Mod%i]*1ll
				*(Mod-Mod/i)%Mod;
		}
		for(int i=2; i<MaxN; ++i)
			inv[i] = 1ll*inv[i-1]*inv[i]%Mod; 
		jc[0] = inv[0] = 1;
	}
	int getC(int n,int m){
		if(m < 0 || n < m) return 0;
		return 1ll*jc[n]*inv[m]%Mod*inv[n-m]%Mod;
	}
	
	int phi[MaxN];
	vector< int > primes;
	void importPhi(){
		phi[1] = 1; // important
		for(int i=2; i<MaxN; ++i)
			phi[i] = i-1;
		for(int i=2,len=0; i<MaxN; ++i){
			if(phi[i] == i-1) ++ len,
				primes.push_back(i);
			for(int j=0; j<len&&primes[j]<=(MaxN-1)/i; ++j){
				phi[i*primes[j]] = phi[i]*(primes[j]-1);
				if(i%primes[j] == 0){
					phi[i*primes[j]] = phi[i]*primes[j];
					break; // till next time
				}
			}
		}
	}
	int getPhi(int n){ return phi[n]; }
}

namespace HardCalc{
	int getPhi(int n){
		return MyMath::getPhi(n);
	}
	// n category, m in total
	inline int getDP(int n,int m){
		if(!n && !m) return 1;
		int res = 0;
		for(int i=0; i<=m/(k+1); ++i){
			int t = MyMath::getC(n,i)*1ll
			*MyMath::getC(m-i*(k+1)+n-1,n-1)%Mod;
			if(i&1) res = (res+Mod-t)%Mod;
			else res = (res+t)%Mod;
		}
		return res;
	}
	int query(int d){
		if((1ll*m*d)%n != 0) return 0;
		int res = 0, t = 1ll*m*d/n;
		for(int i=0; i<=k&&i<=t; ++i)
			res = (res+(i+1ll)*
				getDP(d-t-1,t-i))%Mod;
		return res;
	}
	void solve(){
		int ans = 0, i;
		for(i=1; i<n/i; ++i)
			if(n%i == 0){
				ans = (ans+1ll*getPhi(n/i)*query(i))%Mod;
				ans = (ans+1ll*getPhi(i)*query(n/i))%Mod;
			}
		if(1ll*i*i == n) // sqrt = int
			ans = (ans+1ll*getPhi(i)*query(i))%Mod;
		ans = 1ll*ans*qkpow(n,Mod-2)%Mod;
		printf("%d\n",ans);
	}
}

int main(){
	freopen("gift.in","r",stdin);
	freopen("gift.out","w",stdout);
	MyMath::importC();
	MyMath::importPhi();
	for(int T=readint(); T; --T){
		n = readint(), m = readint();
		k = readint();
		if(!m){ puts("1"); continue; }
		if(m == n){
			putchar('0' + (k == m));
			putchar('\n'); continue;
		}
		HardCalc::solve();
	}
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值