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;
}