轉載自https://www.cnblogs.com/hadilo/p/5914302.html
一、歐幾里得算法(重點是證明,對后續知識有用)
歐幾里得算法,也叫輾轉相除,簡稱 gcd,用於計算兩個整數的最大公約數
定義 gcd(a,b) 為整數 a 與 b 的最大公約數
引理:gcd(a,b)=gcd(b,a%b)
先把這個引理的結論解釋一下:
- 假設a=mc,b=nc,如果a,b的最大公約數是c,則m,n一定互質;
- 假設r=a%b=a-pb=mc-pnc=(m-pn)c,由1可知,r,b的最大公約數是c,且m-pn,n互質,
所以gcd(a,b)==gcd(b,a%b)
證明:
設 r=a%b , c=gcd(a,b)
則 a=xc , b=yc , 其中x , y互質
r=a%b=a-pb=xc-pyc=(x-py)c
而b=yc
可知:y 與 x-py 互質
證明:
假設 y 與 x-py 不互質
設 y=nk , x-py=mk , 且 k>1 (因為互質)
將 y 帶入可得
x-pnk=mk
x=(pn+m)k
則 a=xc=(pn+m)kc , b=yc=nkc
那么此時 a 與 b 的最大公約數為 kc 不為 k
與原命題矛盾,則 y 與 x-py 互質
因為 y 與 x-py 互質,所以 r 與 b 的最大公約數為 c
即 gcd(b,r)=c=gcd(a,b)
得證
當a%b=0時,gcd(a,b)=b
這樣我們可以寫成遞歸形式
int gcd(int a, int b)//最大公約數 { return b == 0 ? a : gcd(b, a % b); } int lcm(int a, int b)//最小公倍數 { return a / gcd(a, b) * b; }
二、擴展歐幾里得算法
擴展歐幾里得算法,簡稱 exgcd,一般用來求解不定方程,求解線性同余方程,求解模的逆元等
引理:存在 x , y 使得 gcd(a,b)=ax+by
證明:
當 b=0 時,gcd(a,b)=a,此時 x=1 , y=0
當 b!=0 時,
設 ax1+by1=gcd(a,b)=gcd(b,a%b)=bx2+(a%b)y2
又因 a%b=a-a/b*b (參見文章開頭的兩個假設,這里的a/b取的是整數(假設里面的p),a-a/b*b!=0)
則 ax1+by1=bx2+(a-a/b*b)y2
ax1+by1=bx2+ay2-a/b*by2
ax1+by1=ay2+bx2-b*a/b*y2
ax1+by1=ay2+b(x2-a/b*y2)
解得 x1=y2 , y1=x2-a/b*y2
因為當 b=0 時存在 x , y 為最后一組解,gcd(a,b)=a,此時 x=1 , y=0
而每一組的解可根據后一組得到
所以第一組的解 x , y 必然存在
得證
根據上面的證明,在實現的時候采用遞歸做法
先遞歸進入下一層,等到到達最后一層即 b=0 時就返回x=1 , y=0
再根據 x=y’ , y=x’-a/b*y’ ( x’ 與 y’ 為下一層的 x 與 y ) 得到當層的解
不斷算出當層的解並返回,最終返回至第一層,得到原解
這里的exgcd求得的x是 ax+by=1的解 (即a,b互質,gcd(a,b)==c)
void exgcd(LL a, LL b, LL &x, LL &y) //拓展歐幾里得算法 { if(!b) x = 1, y = 0; else { exgcd(b, a % b, y, x); y -= x * (a / b); } }
擴展歐幾里德求逆元
LL niYuan(LL a, LL b) //求a對b取模的逆元 { LL x, y; exgcd(a, b, x, y); return (x + b) % b; }
模板題:https://www.cnblogs.com/-citywall123/p/10693036.html
三、exgcd 解不定方程
---------(使用不將a與b轉為互質的方法)-----(gcd(a,b)!=c)
對於 ax+by=c 的不定方程,設 r=gcd(a,b)
當 c%r!=0 時無整數解
當 c%r=0 時,將方程右邊 *r/c 后轉換為 ax+by=r 的形式
可以根據擴展歐幾里得算法求得一組整數解 x0 , y0
而這只是轉換后的方程的解,原方程的一組解應再 *c/r 轉變回去
(如 2x+4y=4 轉換為 2x+4y=2 后應再將解得的 x , y 乘上2)
則原方程解為 x1=x0*c/r , y1=y0*(c/r);
通解 x=x1+b/r*t , y=y1-a/r*t ,其中 t 為整數
證明:
將 x , y 帶入方程得
ax+ab/r*t+by-ab/r*t=c
ax+by=c
此等式恆成立
得證
這里 b/r 與 a/r 為最小的系數,所以求得的解是最多最全面的
證明:
為了推出證明中的 ax+by=c ,且想達到更小的系數,只能將 b/r 與 a/r 同除以一個數 s
而 b/r 與 a/r 互質,且 s 為整數,則 s=1 ,不影響通解
那么 b/r 與 a/r 就為最小的系數
得證
ll x,y,r,s; void exgcd(ll a, ll b, ll &x, ll &y) //拓展歐幾里得算法 { if(!b) x = 1, y = 0; else { exgcd(b, a % b, y, x); y -= x * (a / b); } } ll gcd(ll a,ll b) { return b==0?a:gcd(b,a%b); } void BD(ll a,ll b,ll c,ll r) { exgcd(a,b,x,y); x=x*c/r;//得到原方程的解x和y y=y*c/r; }
模板題:https://www.cnblogs.com/-citywall123/p/10698205.html
模板題:http://www.cnblogs.com/hadilo/p/5917173.html
四、exgcd 解線性同余方程
線性同余方程有解的充分必要條件是當且僅當 c 能夠被 a 與 b的最大公約數整除 c%gcd(a,b)==0
關於 x 的模方程 ax%b=c 的解---------要保證x的系數a為正,如果a小於0,等號兩邊乘以-1 eg:https://www.cnblogs.com/-citywall123/p/10811824.html
方程轉換為 ax+by=c 其中 y 一般為非正整數
(PS怎么轉換:ax%b=c%b
假設一個整數y; 因為ax%b=(ax+by)%b 所以(ax+by)%b=c%b --> ax+by=c)
則問題變為用 exgcd 解不定方程
解得 x1=x0*c/r
通解為 x=x1+b/r*t
設 s=b/r (已證明 b/r 為通解的最小間隔,r=gcd(a,b))
則 x 的最小正整數解為 (x1%s+s)%s
同理可得y的最小正整數解為(y1%ss+ss)%ss;------ss=a/r;
證明:
若 x1>0,則 (x1%s+s)%s=x1%s%s+s%s=x1%s=x1-ts (t∈N)
若 x1<0,因在 C++ 里 a%b=-(-a%b)<0 (a<0 , b>0) 如 -10%4=-2
則 (x1%s+s)%s=(-(-x1%s)+s)%s=(-(ts-x1)+s)%s=ts-x1 (t∈N)
即為 x1 通過加或減上若干個 s 后得到的最小正整數解
得證
亦可偽證 x1<0 的情況:設 x1=-5 , s=2
則 (x1%s+s)%s=(-5%2+2)%2=(-1+2)%2=3%2=1
即為 x1 加上 3 個 s 后的到的最小正整數解
ll x,y,r,s; void exgcd(ll a, ll b, ll &x, ll &y) //拓展歐幾里得算法 { if(!b) x = 1, y = 0; else { exgcd(b, a % b, y, x); y -= x * (a / b); } } ll gcd(ll a,ll b) { return b==0?a:gcd(b,a%b); } void T(ll a,ll b,ll c) { r=gcd(a,b); s=b/r; exgcd(a,b,x,y);//得到x0 x=x*c/r; //得到x1 x=(x%s+s)%s; //得到最小正整數解 }
模板題:https://www.cnblogs.com/-citywall123/p/10694160.html
模板題:https://www.cnblogs.com/-citywall123/p/10680585.html
模板題:http://www.cnblogs.com/hadilo/p/5951091.html