【bzoj1485:】【 [HNOI2009]有趣的数列】模任意数的卡特兰数

（上不了p站我要死了，侵权度娘背锅）

Description
我们称一个长度为2n的数列是有趣的，当且仅当该数列满足以下三个条件：
(1)它是从1到2n共2n个整数的一个排列{ai}；
(2)所有的奇数项满足a1 < a3 < … < a2n-1，所有的偶数项满足a2 < a4 < … < a2n；
(3)任意相邻的两项a2i-1与a2i(1≤i≤n)满足奇数项小于偶数项，即：a2i-1 < a2i。
现在的任务是：对于给定的n，请求出有多少个不同的长度为2n的有趣的数列。因为最后的答案可能很大，所以只要求输出答案 mod P的值。
Input
输入文件只包含用空格隔开的两个整数n和P。输入数据保证，50%的数据满足n≤1000，100%的数据满足n≤1000000且P≤1000000000。
Output
仅含一个整数，表示不同的长度为2n的有趣的数列个数mod P的值。
Sample Input
3 10
Sample Output
5
对应的5个有趣的数列分别为（1,2,3,4,5,6），（1,2,3,5,4,6），（1,3,2,4,5,6），（1,3,2,5,4,6），（1,4,2,5,3,6）。

首先这是一个卡特兰数

我们思考把1~n的数字按顺序摆放。因为分了奇数和偶数，就相当于分了两组队列，我们把这n个数按顺序放入两组队列里。由于任何时候奇数队列里的元素个数都大于等于偶数队列里的，所以是卡特兰的一个经典模型。

答案即为C(2n,n)/(n+1) mod p

但是发现p并不是一个质数，所以不能求逆元，就不能直接用阶乘的公式来做。

这里就有一个很重要的方法了
虽然是模任意数，但是n的范围是o(n)可过的。我们筛出需要的数（如n）的阶乘（n！）中每个素数的指数。具体做法为：
设cnt[i]表示该阶乘（n！）中包含多少个i，进行线性递推。初始值cnt[i]=1。
依次枚举n到1，跳过质数。设minp[i]为i的最小质因子，cnt[i/minp[i]]+=cnt[i], cnt[minp[i]]+=cnt[i], cnt[i]=0。
这部分的代码（init()为线性筛）：

init();
for(int i=n*2;i>=2;i--){
    cnt[i]++;
    if(!notp[i]) continue;
    cnt[i/minp[i]]+=cnt[i];
    cnt[minp[i]]+=cnt[i];
    cnt[i]=0;
}

这样就将除法问题转化为了质数的相减问题，得出C(2n,n)/(n+1)中每个质数的指数，再做快速幂。
时间复杂度：根据某素数定理，n中的素数个数约为 n/ln n 个，快速幂是 log n 的。所以乘起来约为o（n）

完整代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long 
#ifdef WIN32
#define RIN "%I64d"
#else 
#define RIN "%lld"
#endif

const int N=1000000+5;

int cnt1[N*2],cnt2[N*2],cnt3[N*2];
int notp[N*2],cntp=0,prime[N*2],minp[N*2];
int n,p;

void init(){
    notp[1]=1;
    for(int i=2;i<=n*2;i++){
        if(!notp[i])
            prime[++cntp]=i;
        for(int j=1;j<=cntp&&i*prime[j]<=n*2;j++){
            notp[i*prime[j]]=1;
            minp[i*prime[j]]=prime[j];
            if(i%prime[j]==0) break;
        }
    }
}
ll power(int a,int b){
    ll rt=1;
    for(;b;b>>=1,a=(a*a)%p)
        if(b&1) rt=(rt*a)%p;
    return rt;
}
int main(){
    scanf("%d%d",&n,&p);
    init();
    for(int i=n*2;i>=2;i--){
        cnt1[i]++;
        if(!notp[i]) continue;
        cnt1[i/minp[i]]+=cnt1[i];
        cnt1[minp[i]]+=cnt1[i];
        cnt1[i]=0;
    }
    for(int i=n;i>=2;i--){
        cnt2[i]++;
        if(!notp[i]) continue;
        cnt2[i/minp[i]]+=cnt2[i];
        cnt2[minp[i]]+=cnt2[i];
        cnt2[i]=0;
    }
    for(int i=n+1;i>=2;i--){
        cnt3[i]++;
        if(!notp[i]) continue;
        cnt3[i/minp[i]]+=cnt3[i];
        cnt3[minp[i]]+=cnt3[i];
        cnt3[i]=0;
    }
    ll ans=1;
    for(int i=1;i<=cntp;i++){
        int tmp=prime[i];//printf("%d %d %d\n",cnt1[tmp],cnt2[tmp],cnt3[tmp]);
        ans=(ans*power(prime[i],cnt1[tmp]-cnt2[tmp]-cnt3[tmp]))%p;
    }
    printf(RIN"\n",ans);
    return 0;
}