中国剩余定理

MachoMan_1069

于 2019-03-18 13:02:00 发布

阅读量158

点赞数

原文链接：http://www.cnblogs.com/Vscoder/p/10551544.html

版权

前言：

中国剩余定理（$CRT$），也称孙子定理，原文如下：

“有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何?”

很明显这是一个同余方程组，于是我们就可以用中国剩余定理求解

正文：

中国剩余定理

中国剩余定理适用于求解模数两两互质时的同余方程组

设 $b_1,b_2,\ldots,b_n$ 两两互质

则同余方程组 $\begin{cases} x\equiv a_1\ (mod\ b_1) \\ x\equiv a_2\ (mod\ b_2) \\ \qquad\dots \\ x\equiv a_n\ (mod\ b_n) \end{cases}$ 有整数解

并且在 $mod\ lcm=\prod_{i=1}^{n}b_i$ 意义下有唯一解

为 $x=\sum_{i=1}^ka_i*left*inv\;\%\;lcm$

其中 $left=\dfrac{lcm}{b_i}$ ，$inv$ 为 $left\ mod\ b_i$ 意义下的逆元

同时 $CRT$ 一般会配合龟速乘，防止乘起来会爆 $long\ long$

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;
}

ll qmul(ll a,ll b,ll p)
{
    ll ans=0,base=a%p;
    while(b)
    {
        if(b&1) ans=(ans+base)%p;
        base=(base+base)%p;
        b>>=1;
    }
    return ans;
}

ll CRT(int n,ll *a,ll *b)
{
    ll ans=0,lcm=1;
    for(int i=1;i<=n;i++)
        lcm*=b[i];
    for(int i=1;i<=n;i++)
    {
        ll inv,k;
        ll left=lcm/b[i];
        exgcd(left,b[i],inv,k);
        inv=(inv%b[i]+b[i])%b[i];
        ans=(ans+qmul(qmul(a[i],left,lcm),inv,lcm))%lcm;
    }
    return (ans+lcm)%lcm;
}

扩展中国剩余定理

扩展中国剩余定理可以用来解决模数不互质的情况

即同余方程组 $\begin{cases} x\equiv a_1\ (mod\ b_1) \\ x\equiv a_2\ (mod\ b_2) \\ \qquad\dots \\ x\equiv a_n\ (mod\ b_n) \end{cases}$ 中

$b_1,b_2,\ldots,b_n$ 为不一定两两互质的整数，求 $x$ 的最小非负整数解

我们可以假设已经求出了前 $n-1$ 个方程组为一个解为 $x$

且有 $lcm$ 等于前 $n-1$ 个模数的最小公倍数

则前 $n-1$ 个方程组的通解为 $x+k*lcm$

那对与我们加入第 $n$ 个方程后，就是要求一个正整数 $k$

使得 $x+k*lcm \equiv a_n\ (mod\ b_n)$

我们将它稍微整理一下

$lcm*k \equiv a_n-x\ (mod\ b_n)$

设 $a=lcm，b=b_n$

再设 $c=((a[n]-x)\%b+b)\%b$

所以我们就得到了 $ak\equiv c\ (mod\ b)$

也就转化为了解不定方程 $ak+by=c$

我们可以先用 $EXGCD$ 解出 $ak'+by'=gcd(a,b)$

则有$a\;\dfrac{k'}{gcd(a,b)}+b\;\dfrac{y'}{gcd(a,b)}=1$

所以$a\;\dfrac{k'c}{gcd(a,b)}+b\;\dfrac{y'c}{gcd(a,b)}=c$

所以 $k=\dfrac{k'c}{gcd(a,b)}$

所以加入第 $k$ 个同余方程的一个解 $x_n=x+k*lcm$

这样看来，我们可以通过 $k$ 次 $EXGCD$ 来求解这个同余方程组

不会 $EXGCD$ 的可以戳这里

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    ll gcd=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return gcd;
}

ll qmul(ll a,ll b,ll p)
{
    ll ans=0,base=a%p;
    while(b)
    {
        if(b&1) ans=(ans+base)%p;
        base=(base+base)%p;
        b>>=1;
    }
    return ans;
}

ll EXCRT(int n,ll *a,ll *b)
{
    ll ans=a[1],lcm=b[1];
    for(int i=2;i<=n;i++)
    {
        ll k,y;
        ll gcd=exgcd(lcm,b[i],k,y);
        ll c=((a[i]-ans)%b[i]+b[i])%b[i];
        k=qmul(k,c/gcd,b[i]);
        ans+=k*lcm;
        lcm*=b[i]/gcd;
        ans=(ans%lcm+lcm)%lcm;
    }
    return ans;
}