扩展CRT&&扩展lucas

最新推荐文章于 2022-11-02 09:33:05 发布

stone41123

最新推荐文章于 2022-11-02 09:33:05 发布

阅读量600

点赞数

CC 4.0 BY-SA版权

分类专栏：组合数学那本书好厚啊数论真的太难了

本文链接：https://blog.youkuaiyun.com/stone41123/article/details/78877131

数论真的太难了同时被 2 个专栏收录

41 篇文章

订阅专栏

组合数学那本书好厚啊

5 篇文章

订阅专栏

本文提供了扩展中国剩余定理(CRT)与扩展Lucas定理的代码实现。通过这两段代码，可以解决数论中关于同余方程组求解及组合数计算的问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

直接上板子，这里没有证明，只是提供一个还不错的板子
扩展CRT：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstdlib>
#define ll long long
using namespace std;
inline int read(){
    int x=0;char ch=' ';ll f=1;
    while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
    if(ch=='-')f=-1,ch=getchar();
    while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    return x*f;
}
const int N=1005;
inline ll gcd(ll a,ll b){
    if(!b)return a;
    return gcd(b,a%b);
}
inline void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){x=1,y=0;return;}
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
inline ll inv(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    if(x<0)x+=p;
    return x;
}
ll n,c[N],m[N];
inline ll CRT(){
    for(int i=2;i<=n;i++){
        ll m1=m[i-1],m2=m[i],c1=c[i-1],c2=c[i];
        ll t=gcd(m1,m2);
        m[i]=m1*m2/t;
        c[i]=inv(m1/t,m2/t)*((c2-c1)/t)%(m2/t)*m1+c1;
        c[i]=(c[i]%m[i]+m[i])%m[i];
    }
    return c[n];
}
int main(){
    freopen("CRT3.in","r",stdin);
    freopen("CRT3.out","w",stdout);
    n=read();
    for(int i=1;i<=n;i++)c[i]=read();
    for(int i=1;i<=n;i++)m[i]=read();
    printf("%lld",CRT());
    return 0;
}

扩展lucas：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstdlib>
#define ll long long
using namespace std;
inline ll read(){
    ll x=0;char ch=' ';ll f=1;
    while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
    if(ch=='-')f=-1,ch=getchar();
    while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    return x*f;
}
inline void print(ll x){
    if(!x)return;
    print(x/10);
    putchar(x%10+'0');
}
inline void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){x=1,y=0;return;}
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
inline ll inv(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    if(!x)x+=p;
    return x;
}
inline ll ksm(ll a,ll b,ll p){
    ll ans=1;
    while(b){
        if(b&1)ans=(ans*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return ans;
}
inline ll exlucas(ll n,ll pi,ll pk){
    if(!n)return 1LL;
    ll ans=1LL;
    for(ll i=2LL;i<=pk;i++)if(i%pi)ans=(ans*i)%pk;
    ans=ksm(ans,n/pk,pk);
    for(ll i=2LL;i<=n%pk;i++)if(i%pi)ans=(ans*i)%pk;
    return (ans*exlucas(n/pi,pi,pk))%pk;
}
inline ll C(ll n,ll m,ll pi,ll pk){
    if(m>n)return 0LL;
    ll a=exlucas(n,pi,pk),b=exlucas(n-m,pi,pk),c=exlucas(m,pi,pk);
    ll ans,k=0;
    for(ll i=n;i;i/=pi)k+=i/pi;
    for(ll i=n-m;i;i/=pi)k-=i/pi;
    for(ll i=m;i;i/=pi)k-=i/pi;
    ans=a*inv(b,pk)%pk*inv(c,pk)%pk*ksm(pi,k,pk)%pk;
    return ans;
}
ll T,n,m,p,P,ans;
int main(){
    T=read();
    while(T--){
        n=read();m=read();P=read();p=P;ans=0;
        for(ll pi=2;pi*pi<=p;pi++){
            if(p%pi==0){
                ll pk=1LL;
                while(p%pi==0){p/=pi;pk*=pi;}
                ans=(ans+C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P;
            }
        }
        if(p>1)ans=(ans+C(n,m,p,p)*(P/p)%P*inv(P/p,p)%P)%P;
        printf("%lld\n",ans);
    }
    return 0;
}