擴展歐幾里得&&中國剩余定理(學習筆記)


原本是想把CRT、擴展CRT、歐幾里得、擴展歐幾里得都寫在這,但由於博主太菜,剛剛才會EXCRT qwq

在退組邊緣徘徊的我果然還是菜了一點啊!!!

布吉島為什么但就是想奶一口gql省隊穩了2333

不閑扯了,進入正題!


歐幾里得(gcd)&&擴展歐幾里得(exgcd):

先來一個眾人皆知的歐幾里得算法\(gcd\ (\ a,b\ )=gcd\ (\ b,a\ mod\ b\ )\)

證明過程自己完成,這里不多加敘述

再來普及一個簡單的裴蜀定理:若a,b是整數,且 \(gcd(a,b)=d\),那么一定存在整數 \(x,y\),使得\(ax+by=d\) 成立。

裴蜀定理證明請看這里(主要是因為我之前打了一段后發現沒證完然后就懶得打了qwq)

那么 \(\bf\Huge\text{擴展歐幾里得}​\) 來了!!!

擴展歐幾里得算法:可以用來求解 \(ax+by=gcd(a,b)\) ,也就是求同余方程 \(ax\equiv gcd(a,b)(mod\ b)\) 的解

算法過程及解如下:

我們要求解 \(ax+by=gcd(a,b)\) 其實主要是想算出 \(x\) 的值,\(y\) 是一個輔助解

然后我們構造這么一個式子:\(bx_1+(a\ mod\ b)y_1=gcd(b,(a\ mod\ b))\) ,根據歐幾里得算法,我們可得出 \(bx_1+(a\ mod\ b)y_1=ax+by\) ,不妨把 \(a\ mod\ b\) 變成 \(a-b*(a/b)\) ,然后帶入原式:


\(bx_1+(a-b*(a/b))y_1=ax+by\)

\(bx_1+ay_1-b*(a/b)y_1=ax+by\)

\(ay_1+b(x_1-(a/b)y1)=ax+by\)


那么我們求出了一組解:\(x=y_1,y=x_1-(a/b)y_2\)

然后將我們求解的兩個式子對比一下:


\(ax+by=gcd(a,b)\)

\(bx_1+(a\ mod\ b)y_1=gcd(b,(a\ mod\ b))\)


是不是發現了什么?是不是有點像歐幾里得算法?

然后我們按這種方式遞推下去,直到 \(b=0\) ,那么:


\(ax+by=gcd(a,b)\)

\(ax=gcd(a,0)\)


顯然 \(gcd(a,0)=a,x=1\) ,此時已經和 \(y\) 的值沒有關系了,即此時 \(y\) 的值是任意數,但是這里建議把 \(y\) 賦成 \(0\) ,可以避免在返回中爆 \(long\ long\)

之后我們已經求出來了關於同余方程的一個解 \(x\) ,雖然 \(x\) 不一定是最小的,但顯然 \(x\) 加上或減去 \(b\) 是沒有任何影響的,所以用一個 \(x = (x \% b + b) \% b\) 就可以了求出滿足同余方程的最小正整數解了

上代碼:

void exgcd(int a,int b,int &x,int &y)
{//x,y都是要返回的值
    if(b==0)
    {
        x=1,y=0;
        return;
    }
    exgcd(b,(a%b),x,y);//遞歸
    int xx=y;//記錄一下解
    y=x-(a/b)*xx,x=xx;
    return;
}

至此,擴展歐幾里得講解完畢qwq


中國剩余定理(CRT):

中國剩余定理其實真的不難qwq

一般情況下我們是要求解如下式子:

\(\begin{cases} x\equiv a_1(\mod b_1)\quad \\ x\equiv a_2(\mod b_2)\quad \\ ...\quad \\ x\equiv a_n(\mod b_n)\quad \\ \end{cases}\)

其中所有的 \(b_1,b_2...b_n​\) 都互質,這里的 \(a_i​\) 都小於 \(b_i​\)

首先我們構造一個數 \(B=b_1*b_2*...*b_n​\) ,那么顯然當我們求出 \(x​\) 之后加上或減去 \(B​\) 都是成立的。

接下來考慮對於每一個同余方程的處理:

對於同余式 \(x\equiv a_k(\mod b_k)\) ,我們可以構造一個 \(B_k=B/b_k\) ,顯然 \(gcd(B_k,b_k)=1\) ,根據裴蜀定理可得:存在整數 \(i,j\) 使得 \(iB_k+jb_k=1\) ,即 \(iB_k\equiv 1(\mod b_k)\) 。因為 \(iB_k\ mod\ b_k=1\) ,那么有 \(a_i*iB_k\ mod\ b_k=a_i\)\(i\) 可以通過擴展歐幾里得求得,所以該同余方程的一個解是 \(x=a_i*i*B_k\)

再回到整個方程組,我們繼續構造一個數 \(x=a_1i_1B_1+a_2i_2B_2+...+a_ni_nB_n\) ,那么這個數就是同余方程組的一個解。因為對於第 \(k\) 個同余方程,\(B_{1...n}\) 中除 \(B_k\) 之外所有數都是 \(b_k\) 的倍數,所以同余方程是成立的。

放一下CRT的代碼:

void exgcd(int a,int b,int &x,int &y){
    if(b==0)
    {
        x=1,y=0;
        return;
    }
    exgcd(b,(a%b),x,y);
    int xx=y;
    y=x-(a/b)*xx,x=xx;
    return;
}//擴展歐幾里得求解同余方程

lt china()
{
    for(int i=1;i<=n;i++) N*=B[i];
    for(int i=1;i<=n;i++)
    {
        int a=N/B[i],b=B[i],x=0,y=0;
        exgcd(a,b,x,y);//此處的x還要變成最小正整數解
        ans+=a*((x%b+b)%b)*A[i];
    }
    return (ans%N+N)%N;
}

擴展中國剩余定理(EXCRT)

好了,同樣是形如下面的式子:

\(\begin{cases} x\equiv a_1(\mod b_1)\quad \\ x\equiv a_2(\mod b_2)\quad \\ ...\quad \\ x\equiv a_n(\mod b_n)\quad \\ \end{cases}​\)

只不過這次不保證所有的 \(b_1,b_2,...,b_n\) 互質

因為不互質,所以我們無法像CRT一樣使用擴歐。

考慮你已經求出來了前 \(k-1​\) 個同余方程的解,得到的解為 \(ans​\) ,設 \(lcm​\) 為前 \(k-1​\) 個方程中所有的 \(b_i​\) 的最小公倍數,則前 \(k-1​\) 個方程的解為 \(x=ans+i*lcm​\) ,而我們只需要確定一個 \(i​\) 使得 \(ans+i*lcm\equiv a_k(\mod b_k)​\) ,然后更新一下 \(lcm​\) 就好了。

又是一波轉化:


\(ans+i*lcm\equiv a_k(\mod b_k)\)

\(i*lcm\equiv a_k-ans(\mod b_k)\)


注意:此處的 \(lcm\)\(b_k\) 不一定互質,需要稍加處理,設 \(gcd=gcd(lcm,b_k),c=(a_k-ans)\ mod\ b_k\) ,由上式可得:

\(i*lcm+h*b_k=c\)


由裴蜀定理可得此方程有解的必要條件是 \(gcd|c\) ,於是繼續變形:


\(i/gcd*lcm+h*b_k/gcd=c/gcd\)

\(i/gcd*lcm\equiv c/gcd(\mod b_k/gcd)\)


此時 \(gcd(lcm,b_k/gcd)=1\) ,對於上面這個方程,我們依然可以先求出一個 \(j\) 滿足 \(j*lcm\equiv 1(\mod b_k/gcd)\) ,然后再乘以 \(c/gcd\) 倍就好了。

具體代碼中有解釋:

int exgcd(int a,int b,int &x,int &y)
{//擴展歐幾里得,一並求出gcd
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    int res=exgcd(b,(a%b),x,y);
    int xx=x;
    x=y,y=xx-(a/b)*y;
    return res;
}

int mul(int s,int p,int mod)
{//這是一個快速乘,防止s*p爆long long
    int res=0;
    while(p)
    {
        if(p%2) res=(res+s)%mod;
        s=(s+s)%mod;
        p/=2;
    }
    return res%mod;
}

int EXCRT()
{
    m=b[1],ans=a[1];//第一個方程要特殊處理,直接賦值就好了
    for(int i=2;i<=n;i++)
    {
        int x,y,c=(a[i]-ans%b[i]+b[i])%b[i];
        //c就是前面的ak-ans
        int gcd=exgcd(m,b[i],x,y);
        x=mul(x,c/gcd,b[i]/gcd);
        //這一段我也弄了好久,最后終於搞懂了
        //大致把x為什么要乘c/gcd,和%(b/gcd)的原因寫在上面了
        //不懂歡迎提問
        ans+=x*m;//更新ans
        m*=(b[i]/gcd);//更新lcm
        ans=(ans%m+m)%m;
    }
    return (ans%m+m)%m);
}

終於更完了qwq心累


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM