1. 程式人生 > 實用技巧 >題解 P3200 【[HNOI2009]有趣的數列】

題解 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=\binom{2n}{n}-\binom{2n}{n-1}\)

問題變為求 \(\binom{n}{m}\;mod\;p\),但也許是出於某些不可告人的因素,這個 \(p\) 不是一個質數。

考慮擴充套件 \(\text{Lucas}\) 定理之後用 \(\text{CRT}\) 合併,但是我們需要把 \(p\) 分解成若干個質數的冪的積。所以如果 \(p\)

是一個質數的話,按照擴充套件 \(\text{Lucas}\) 的方法複雜度是 \(O(p)\) 的。

但是 \(n\) 的範圍不大,我們可以根號分治。對於 \(p\),它最多隻可能有一個大於 \(\sqrt{p}\) 的質因子。

設這個質因子為 \(x\),如果 \(x>n\),那麼在模 \(x\) 意義下,\(1-n\) 這些數的階乘都是有逆元的,直接算就可以了。如果 \(x \leq n\),我們再擴充套件 \(\text{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\n",catalan(n));	
	return 0;
}