题解 P3200 【[HNOI2009]有趣的数列】

题目链接

Solution [HNOI2009]有趣的数列

题目大意：求有多少个 (1-2n) 的排列，满足(a_1 < a_3 < a_5 < cdots < a_{2n-1})，(a_2 < a_4 < a_6 < cdots < a_{2n}),且 (forall k,a_{2k-1}<a_{2k})

扩展卢卡斯

---

分析：

生成这个排列，相当于在数 (1-2n) 中选出 (n) 个数，顺次放入奇数项，剩下的 (n) 个数再顺次放入偶数项。

这样满足了前两个要求，对于第三个要求，我们把 (forall a_{2k-1},a_{2k}) 称为一对。

那么排列满足第三个要求，就相当于顺次从 (1) 枚举到 (2n)，如果这个数要丢进偶数项，我们就将它丢到最后一个没有配对的奇数项的后面，且这 (n) 步操作都合法。

然后这其实就是求合法括号序列，于是答案为 (Cat_n=inom{2n}{n}-inom{2n}{n-1})

问题变为求 (inom{n}{m};mod;p)，但也许是出于某些不可告人的因素，这个 (p) 不是一个质数。

考虑扩展 ( ext{Lucas}) 定理之后用 ( ext{CRT}) 合并，但是我们需要把 (p) 分解成若干个质数的幂的积。所以如果 (p) 是一个质数的话，按照扩展 ( ext{Lucas}) 的方法复杂度是 (O(p)) 的。

但是 (n) 的范围不大，我们可以根号分治。对于 (p)，它最多只可能有一个大于 (sqrt{p}) 的质因子。

设这个质因子为 (x)，如果 (x>n)，那么在模 (x) 意义下，(1-n) 这些数的阶乘都是有逆元的，直接算就可以了。如果 (x leq n)，我们再扩展 ( ext{Lucas})。

分解 (p) 的时候同样要根号分治，可能只有我犯这么憨憨的错吧

#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
int exgcd(const int a,const int b,ll &x,ll &y){
	if(!b){
		x = 1,y = 0;
		return a;
	}
	int res = exgcd(b,a % b,y,x);
	y -= (a / b) * x;
	return res;
}
int inv(const int a,const int mod){
	ll x,y;
	exgcd(a,mod,x,y);
	return ((x % mod) + mod) % mod;
}
int n,mod;
namespace lucas{
	int mod,p;
	int add(const int a,const int b){return (a + b) % mod;}
	int mul(const int a,const int b){return (1ll * a * b) % mod;}
	int qpow(int base,int b){
		int res = 1;
		while(b){
			if(b & 1)res = mul(res,base);
			base = mul(base,base);
			b >>= 1;
		}
		return res;
	}
	int f(const int n){
		if(n == 0)return 1;
		int rou = 1,rem = 1;
		for(int i = 1;i <= mod;i++)
			if(i % p)rou = mul(rou,i);
		rou = qpow(rou,n / mod);
		for(ll i = mod * (n / mod) + 1;i <= n;i++)
			if(i % p)rem = mul(rem,i % mod);
		return mul(f(n / p),mul(rou,rem));
	}
	int calc(const int n){
		ll res = 0;
		for(int now = p;now <= n;now *= p)
			res += (n / now);
		return res;
	}
	int C(const int n,const int m){
		int res = f(n);
		res = mul(res,inv(mul(f(m),f(n - m)),mod));
		res = mul(res,qpow(p,calc(n) - calc(m) - calc(n - m)));
		return res;
	}
};
namespace crt{
	const int maxn = 32;
	int a[maxn],m[maxn],tot;
	int crt(){
		int M = 1;
		for(int i = 1;i <= tot;i++)M = M * m[i];
		int res = 0;
		for(int i = 1;i <= tot;i++){
			const int tmp = M / m[i];
			const int inv = ::inv(tmp,m[i]);
			res += ((ll)a[i] * ((ll)tmp * inv % M) % M);
			res %= M;
		}
		return res;
	}
}
int fac(const int n,const int p){
	int res = 1;
	for(int i = 1;i <= n;i++)res = ((ll)res * i) % p;
	return res;
}
int C(const int n,const int m){
	int p = mod;
	crt::tot = 0;
	const int lim = sqrt(n) + 10;
	for(int i = 2;i <= lim;i++)
		if(p % i == 0){
			lucas::p = i;
			lucas::mod = 1;
			do{
				lucas::mod *= i;
				p /= i;
			}while(p % i == 0);
			crt::tot++;
			crt::m[crt::tot] = lucas::mod;
			crt::a[crt::tot] = lucas::C(n,m);
		}
	if(p != 1){
		lucas::p = lucas::mod = p;
		crt::tot++;
		crt::m[crt::tot] = p;
		if(p > n)crt::a[crt::tot] = ((ll)fac(n,p) * inv((ll)fac(m,p) * fac(n - m,p) % p,p)) % p;
		else crt::a[crt::tot] = lucas::C(n,m);
	}
	return crt::crt();
}
int catalan(const int n){
	return (C(n << 1,n) - C(n << 1,n - 1) + mod) % mod;
}
int main(){
#ifndef ONLINE_JUDGE
	freopen("fafa.in","r",stdin);
#endif
	scanf("%d %d",&n,&mod);
	printf("%d
",catalan(n));	
	return 0;
}