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

最新推荐文章于 2024-12-04 22:23:47 发布

原创最新推荐文章于 2024-12-04 22:23:47 发布 · 369 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#CRT #中国剩余定理

本文探讨了在非质数模意义下杨辉三角形的计算问题，利用中国剩余定理(CRT)进行求解，并详细介绍了算法实现过程，包括模幂运算、逆元计算等关键技术。

这里写图片描述

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