NOIP模拟34

考试的时候被T2卡了一年。。。。考虑了一下正解的式子，然后没去给左边分解因数，去给后面乘倍数。。。理论复杂度O（n^2），实际好像卡不掉的样子。但是由于我智障的打了一棵主席树，他M了。。。。

预计得分100+？？+20，实际得分100+70+20

T3，

　　这道题dp式子想一想就出来了，但是由于模数不保证质数，如果用组合数要exlucas，其实可以不用组合数，但是由于我过于智障，还是打了组合数打法，但是由于我不会懒得打ex，于是我们可以。。。。。分解质因数。。。。。

　　首先求出1到m每个数的质因子，然后对于组合数分解成质因数相乘的形式，如果每求一个组合数都暴力将所有质因子乘一遍的话，显然会死，考虑优化：

　　从当前组合数递推到下一个质组合数只需要乘一个数，再除掉一个数，我们观察到除掉的那个数必然<5000，于是我们可以每次将小于5000的质因子暴力乘一遍，再乘上大于5000的质因子的乘积，那么复杂度就对了。

　　注意在求每个数的质因子时，我们不需要每个数O(sqrt(n))求，可以在线筛的同时处理出来，具体来说，i×pri的质因子一定只比i多一个pri，我们可以记录每个数的前趋是谁，这样复杂度O（n），也可以直接让乘积继承i的vector，复杂度和埃氏筛一样，为O（nloglogn），都是可以接受的。

代码比exlucas短，效率也快很多（虽然不如不用组合数）

#include<bits/stdc++.h>
#define ll long long
using namespace std;
int n,m,p,pri[1000005],tot,mx,a[1000005];
vector<int>h[1000005];
char v[1000005];
ll num[5005][5005],f[2][5005],c[5005],aa=1,sum[1000005];
inline int read(int x=0,char ch=getchar()){
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9')x=x*10-48+ch,ch=getchar();
    return x;
}
inline ll qpow(ll x,ll y,ll ans=1){
    for(;y;y>>=1,x=x*x%p)
        if(y&1)
            ans=ans*x%p;
    return ans;
}
signed main(){
    scanf("%d%d%d",&n,&m,&p);
    for(int i=1;i<=n;i++) a[i]=read(),mx=max(mx,a[i]);
    for(int i=2;i<=m;i++){
        if(!v[i]) pri[++tot]=i,h[i].push_back(i);
        for(int j=1;j<=tot&&i*pri[j]<=m;j++){
            v[i*pri[j]]=true;
            h[i*pri[j]]=h[i];
            h[i*pri[j]].push_back(pri[j]);
            if(i%pri[j]==0) break;
        }
    }
    for(int j=0;j<h[m].size();j++){
        if(h[m][j]>5000)aa=aa*h[m][j]%p;
        else sum[h[m][j]]++;
    }
    for(int i=1;i<=min(mx,m);i++){
        c[i]=aa;
        for(int j=1;j<=tot&&pri[j]<=5000;j++)
            c[i]=c[i]*qpow(pri[j],sum[pri[j]])%p;
        if(i!=min(mx,m)){
            for(int j=0;j<h[i+1].size();j++)
                sum[h[i+1][j]]--;
            for(int j=0;j<h[m-i].size();j++)
                if(h[m-i][j]>5000)aa=aa*h[m-i][j]%p;
                else sum[h[m-i][j]]++;
        }
    }
    num[0][0]=1;
    for(int i=1;i<=mx;i++)
        for(int j=1;j<=min(mx,m);j++)
            num[i][j]=(num[i-1][j-1]*j+num[i-1][j]*(j-1))%p;
    for(int i=1;i<=min(a[1],m);i++)
        f[1][i]=num[a[1]][i]*c[i]%p,sum[1]=(sum[1]+f[1][i])%p;
    for(int i=2;i<=n;i++){
        sum[i]=0;
        for(int j=1;j<=min(a[i],m);j++){
            f[i&1][j]=((sum[i-1]*num[a[i]][j]%p*c[j]-(j<=a[i-1])*f[(i-1)&1][j]*num[a[i]][j])%p+p)%p;
            sum[i]=(sum[i]+f[i&1][j])%p;
        }
    }
    printf("%lld\n",sum[n]);
    return 0;
}