[CRT][中国剩余定理]膜法

来源:互联网 发布:类似sketch windows 编辑:程序博客网 时间:2024/05/21 19:25

这里写图片描述

容易发现是杨辉三角形,第i个数的贡献为c(n-1,i-1);但是由于mod数不是质数,所以考虑用CRT来做。先拆mod数,然后在求组合数的过程中将每个数拆成a*pi^ci的形式,非常感谢dmsdalao对我中国剩余定理的指导。
关于CRT详见
http://blog.csdn.net/qq_36993218/article/details/60956475


#include<cstdio>#include<algorithm>#include<iostream>#define ll long long#define Pair pair<ll,ll>using namespace std;int n,K;inline int read(){    char c;    bool flag=false;    while((c=getchar())>'9'||c<'0')    if(c=='-') flag=true;    int res=c-'0';    while((c=getchar())>='0'&&c<='9')    res=(res<<3)+(res<<1)+c-'0';    return flag?-res:res;}const int N=4*1e5+7;int tot,a[N];ll p[N],c[N],m[N],M[N],H[N],ans[N];inline ll ksm(ll a,int b){    ll res=1;    while(b)    {        if(b&1) res=res*a;        a=a*a;        b>>=1;    }    return res;}inline ll ksm(ll a,int b,int pyz){    ll res=1;    while(b)    {        if(b&1) res=res*a%pyz;        a=a*a%pyz;        b>>=1;    }    return res;}inline ll inv(ll a,int x){    return ksm(a,(p[x]-1)*m[x]/p[x]-1,m[x]);}Pair fac[N];inline Pair div(Pair A,Pair B,int x){    return Pair(A.first*inv(B.first,x)%m[x],A.second-B.second);}inline Pair mul(Pair A,Pair B,int x){    return Pair(A.first*B.first%m[x],A.second+B.second);}inline Pair C(int n,int m,int x){    return div(fac[n],mul(fac[m],fac[n-m],x),x);}inline Pair trans(int i,int x){    Pair res=Pair(0,0);    while(i%p[x]==0&&i)    {        i/=p[x];        res.second++;    }     res.first=i;    return res;}inline ll get_num(Pair x,int i){    return (ll)x.first*ksm(p[i],x.second,m[i])%m[i];} inline void solve(int now){    fac[0]=Pair(1,0);    for(int i=1;i<=n;++i)        fac[i]=mul(fac[i-1],trans(i,now),now);    for(int i=1;i<=n;++i)        ans[now]=(ans[now]+get_num(mul(trans(a[i],now),C(n-1,i-1,now),now),now))%m[now];} int main(){    freopen("magic.in","r",stdin);    freopen("magic.out","w",stdout);    n=read();    K=read();    for(int i=1;i<=n;++i) a[i]=read();    int tmp=K;    for(int i=2;i<=tmp;++i)    if(tmp%i==0)    {        p[++tot]=i;        while(tmp%i==0)        {            c[tot]++;            tmp/=i;        }    }    for(int i=1;i<=tot;++i)    {        m[i]=ksm(p[i],c[i]);        M[i]=K/m[i];        H[i]=inv(M[i],i);        solve(i);    }    ll Ans=0;    for(int i=1;i<=tot;++i)        Ans=(Ans+(ll)ans[i]*M[i]%K*H[i]%K)%K;    cout<<Ans; }
0 0