扩展中国剩余定理(EXCRT)

 

数论基础知识吧

主要还是留个板子,之后就不用动脑子了(

模板:

//在主函数依次调用
//fold(a,b,m,n)和excrt(b,m,n)

//求解ax+by=gcd(a,b),返回gcd(a,b) 
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0LL)
    {
        x=1,y=0;
        return a;
    }
    
    ll tmp=exgcd(b,a%b,y,x);
    y=y-a/b*x;
    return tmp;
}

inline ll mod(ll x,ll m)
{
    if(x>=m)
        x-=m;
    if(x<0)
        x+=m;
    return x;
}

//防止中间结果爆longlong的x*y%m 
ll mul(ll x,ll y,ll m)
{
    ll res=0,tmp=x;
    x=mod(x%m+m,m),y=mod(y%m+m,m);
    while(y)
    {
        if(y&1LL)
            res=mod(res+tmp,m);
        tmp=mod(tmp*2,m);
        y/=2;
    }
    return res;
}

//求解方程组x_i=a_i mod m_i 
ll excrt(ll *a,ll *m,int n)
{
    ll res=0,lcm=1;
    for(int i=1;i<=n;i++)
    {
        ll A,B,C,x,y,gcd;
        A=lcm,B=m[i],C=a[i]-res;
        gcd=exgcd(A,B,x,y);
        
        //无解 
        if(C%gcd!=0)
            return -1;
        
        ll tmp=lcm;
        lcm=lcm/gcd*m[i];
        res=res+mul(C/gcd,mul(tmp,x,lcm),lcm);
        res=mod(res%lcm+lcm,lcm);
    }
    return res;
}

//将a_i的逆元乘到b_i上
//之后要excrt的是b_i与m_i 
ll fold(ll *a,ll *b,ll *m,int n)
{
    for(int i=1;i<=n;i++)
    {
        ll x,y,gcd;
        gcd=exgcd(a[i],m[i],x,y);
        
        //无逆元 
        if(b[i]%gcd!=0)
            return -1;
        
        a[i]/=gcd,b[i]/=gcd,m[i]/=gcd;
        exgcd(a[i],m[i],x,y);
        b[i]=mul(b[i],x,m[i]); 
    }
    return 1;
}
View Code

 


 

主要用来求解线性同余方程组

即我们希望对于以下方程组求最小的$x$,使得

\[\begin{equation} \begin{cases} x\equiv a_1 (mod\ \ m_1)\\ x\equiv a_2(mod\ \ m_2)\\ \vdots \\ x\equiv a_n(mod\ \ m_n)\end{cases} \end{equation}\nonumber \]

其中$m_1,m_2,...,m_n$不保证互质

 


 

在此之前需要先掌握扩展欧几里德(exgcd)

参考这篇就足够了,很详细(最后一步式子错了,不过结论对的):RiverHamster - exgcd

exgcd解决的问题是:$ax+by=gcd(a,b)$

//求解ax+by=gcd(a,b),返回gcd(a,b) 
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0LL)
    {
        x=1,y=0;
        return a;
    }
    
    ll tmp=exgcd(b,a%b,y,x);
    y=y-a/b*x;
    return tmp;
}
View Code

 

一些除了本身之外的用法:

1. 解$ax+by=c$的不定方程:将解出的$x,y$同乘$\frac{c}{gcd(a,b)}$即可得到特解;若$c\ mod\ gcd(a,b)\neq 0$,则方程无解

2. 求逆元:若$a,b$互质,那么求出的$x$就是$a$在模$b$意义下的逆元了

3. 解$Ax\equiv B(mod\ \ M)$的同余方程:令$a=A,b=M,c=B$,则$X=x\cdot \frac{c}{gcd(a,b)}$为该方程的解;若$c\ mod\ gcd(a,b)\neq 0$,则该方程无解

 


 

有了上面的铺垫,现在可以愉快地求同余方程组了

参考这篇:niiick - 题解 P4777 【【模板】扩展中国剩余定理(EXCRT)】

假设现在已经算完了前$k-1$个方程组的解,设$M=lcm(m_1,m_2,...,m_{k-1})$

则通解为$x+t\cdot M$($t$为任意整数)

对于第$k$个方程,我们希望$x+t\cdot M\equiv a_k(mod\ \ m_k)$

即,希望找到合适的$t,s$使得$x+t\cdot M=a_k+s\cdot m_k$

移项后就是$t\cdot M+s\cdot m_k=a_k-x$,已经是exgcd可以处理的$ax+by=c$的形式了;用$exgcd(M,m_k,t,s)$,再乘以$\frac{a_k-x}{gcd(M,m_k)}$,即可得到$t,s$

于是前$k$个方程的特解为$x'=x+t\cdot M$,如令$M'=lcm(M,m_k)$则通解为$x'+t'M'$($t'$为任意整数)

这样一步一步做下去,直到算出前$n$个方程的特解

于是可以写出下面的代码(题目为:Luogu P4777 (【模板】扩展中国剩余定理(EXCRT))

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
const int N=100005;

//求解ax+by=gcd(a,b),返回gcd(a,b) 
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0LL)
    {
        x=1,y=0;
        return a;
    }
    
    ll tmp=exgcd(b,a%b,y,x);
    y=y-a/b*x;
    return tmp;
}

inline ll mod(ll x,ll m)
{
    if(x>=m)
        x-=m;
    if(x<0)
        x+=m;
    return x;
}

//防止中间结果爆longlong的x*y%m 
ll mul(ll x,ll y,ll m)
{
    ll res=0,tmp=x;
    x=mod(x%m+m,m),y=mod(y%m+m,m);
    while(y)
    {
        if(y&1LL)
            res=mod(res+tmp,m);
        tmp=mod(tmp*2,m);
        y/=2;
    }
    return res;
}

//求解方程组x_i=a_i mod m_i 
ll excrt(ll *a,ll *m,int n)
{
    ll res=0,lcm=1;
    for(int i=1;i<=n;i++)
    {
        ll A,B,C,x,y,gcd;
        A=lcm,B=m[i],C=a[i]-res;
        gcd=exgcd(A,B,x,y);
        
        //无解 
        if(C%gcd!=0)
            return -1;
        
        ll tmp=lcm;
        lcm=lcm/gcd*m[i];
        res=res+mul(C/gcd,mul(tmp,x,lcm),lcm);
        res=mod(res%lcm+lcm,lcm);
    }
    return res;
}

int n;
ll a[N],m[N];

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&m[i],&a[i]);
    
    printf("%lld\n",excrt(a,m,n));
    return 0;
}
View Code

 

也可以有再稍微扩展一点的做法,即解决

\[\begin{equation} \begin{cases} a_1x\equiv b_1 (mod\ \ m_1)\\ a_2x\equiv b_2(mod\ \ m_2)\\ \vdots \\ a_nx\equiv b_n(mod\ \ m_n)\end{cases} \end{equation}\nonumber \]

当$gcd(a_i,m_i)=1$时,可以通过exgcd求逆元

当$gcd(a_i,m_i)\neq 1$,若$b_i\ mod\ gcd(a_i,m_i)\neq 0$,则该方程必然无解;否则将$a_i,b_i,m_i$同除$gcd(a_i,m_i)$,则$\frac{a_i}{gcd(a_i,m_i)},\frac{m_i}{gcd(a_i,m_i)}$互质,那么也可以通过exgcd求逆元

之前看错题意时想的,没找到题目

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
const int N=100005;

//求解ax+by=gcd(a,b),返回gcd(a,b) 
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0LL)
    {
        x=1,y=0;
        return a;
    }
    
    ll tmp=exgcd(b,a%b,y,x);
    y=y-a/b*x;
    return tmp;
}

inline ll mod(ll x,ll m)
{
    if(x>=m)
        x-=m;
    if(x<0)
        x+=m;
    return x;
}

//防止中间结果爆longlong的x*y%m 
ll mul(ll x,ll y,ll m)
{
    ll res=0,tmp=x;
    x=mod(x%m+m,m),y=mod(y%m+m,m);
    while(y)
    {
        if(y&1LL)
            res=mod(res+tmp,m);
        tmp=mod(tmp*2,m);
        y/=2;
    }
    return res;
}

//求解方程组x_i=a_i mod m_i 
ll excrt(ll *a,ll *m,int n)
{
    ll res=0,lcm=1;
    for(int i=1;i<=n;i++)
    {
        ll A,B,C,x,y,gcd;
        A=lcm,B=m[i],C=a[i]-res;
        gcd=exgcd(A,B,x,y);
        
        //无解 
        if(C%gcd!=0)
            return -1;
        
        ll tmp=lcm;
        lcm=lcm/gcd*m[i];
        res=res+mul(C/gcd,mul(tmp,x,lcm),lcm);
        res=mod(res%lcm+lcm,lcm);
    }
    return res;
}

//将a_i的逆元乘到b_i上
//之后要excrt的是b_i与m_i 
ll fold(ll *a,ll *b,ll *m,int n)
{
    for(int i=1;i<=n;i++)
    {
        ll x,y,gcd;
        gcd=exgcd(a[i],m[i],x,y);
        
        //无逆元 
        if(b[i]%gcd!=0)
            return -1;
        
        a[i]/=gcd,b[i]/=gcd,m[i]/=gcd;
        exgcd(a[i],m[i],x,y);
        b[i]=mul(b[i],x,m[i]); 
    }
    return 1;
}

int n;
ll a[N],b[N],m[N];

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&m[i],&b[i]),a[i]=1;
    
    fold(a,b,m,n);
    printf("%lld\n",excrt(b,m,n));
    return 0;
}
View Code

 


 

只是留个模板,毕竟玩不出什么花里胡哨的

 

(完)

posted @ 2020-01-24 13:13  LiuRunky  阅读(332)  评论(0编辑  收藏  举报