Algorithm: GCD、EXGCD、Inverse Element
數論基礎
數論是純數學的一個研究分支,主要研究整數的性質。初等數論包括整除理論、同余理論、連分數理論。這一篇主要記錄的是同余相關的基礎知識。
取模
取模是一種運算,本質就是帶余除法,運算結果就是餘數。取模運算結果的符號由被模數(被除數)決定。
\[ 7\%4=3;\space7\%(-4)=3;\\ (-7)\%4=-3;\space(-7)\%(-4)=-3 \]
取模運算的性質
\[ 設a>b>0,有:\\ (a+b)\%c=(a\%c+b\%c)\%c\\ (a-b)\%c=(a\%c-b\%c+c)\%c \\ (a\times b)\%c=(a\%c\times b\%c)\%c \]
例題1
\[ 給定兩個非負整數a,b和正整數n,計算f(a^b)\%n。\\ 其中,f(0)=f(1)=1,且對於\forall i>0,f(i+2)=f(i+1)+f(i)\\ (0\le a,b\le 2^{64},1\le n\le 1000) \]
找規律,對於mod n來說,最多n²項就會出現重複,找出重複項即可得到循環周期,然後找對應的項即可。
例題2:hdu 1212
題意是給一個位數不超過1,000的正整數A和一個大小不超過100,000的正整數B,求A%B。
這裏只需要用到上面的兩條性質,把A截斷,從高到低模擬做除法即可。
#include <iostream>
string s; int mod, ans;
int main() {
while(std::cin >> str >> mod) {
ans = 0;
for(int i = 0; s[i]; i++) ans = (ans * 10 + (s[i] - '0')) % mod;
std::cout << ans << std::endl;
}
return 0;
}
GCD
GCD即Greatest Common Divisor,最大公約數。最大公約數即能同時整除給定兩個整數的最大正整數。特殊地:gcd(a,0)=a。求解最大公約數通常使用輾轉相除法。下面是輾轉相除法的代碼實現:
typedef long long LL;
LL gcd(LL a, LL b) {
return b ? gcd(b, a % b) : a;
}
輾轉相除法得到的餘數序列增長速度比斐波那契數列更快,已知斐波那契的增長是指數級別,則輾轉相除法的複雜度是對數級別。
輾轉相除法的證明:
1.設兩個數a、b(a>b),他們的最大公約數為gcd(a,b),r = a % b,k = (a – a % b)/ b,那麼證明輾轉相除法,即是要證明:gcd(a,b)=gcd(b,r)。首先我們令c = gcd(a,b),那麼肯定存在互質的整數m,n,使得a = mc,b = nc。那麼有r = a – kb = mc – knc = (m – kn)c,根據這條式子,c也是r的因數。回過頭看,如果能證明m-kn和n互質,那麼就可以說明c=gcd((m-kn)c,nc)=gcd(b,r),所以把問題再次轉化為:求證m-kn和n互質。
2.反證法證明m-kn和n互質:假設m-kn和n不互質,用數學語言描述為:假設存在整數x,y,d,其中d>1,使得m-kn=xd,n=yd。那麼有m=kn+xd=kyd+xd=(ky+x)d,從而a=mc=(ky+x)cd,b=nc=ycd,則gcd(a,b)=cd≠c,與前設矛盾。故m-kn和n互質得證,也就證明了gcd(a,b)=gcd(b,r)。
唯一分解定理
對於任意大於2的正整數X,它總能寫成如下形式:
\[ X=p_1^{a_1}\times p_2^{a_2} \times ……\times p_k^{a_k},其中p是各不相同的質數 \]
即質因數分解。且這個分解唯一。
LCM
LCM即Least Common Multiple,最小公倍數。
\[ \begin{align} &根據唯一分解定理,對於給定的a,b:\\ &a=p_1^{a_1}\times p_2^{a_2} \times……\times p_k^{a_k}\\ &b=p_1^{b_1}\times p_2^{b_2} \times……\times p_k^{b_k}\\ &那麼:\\ &gcd(a,b)=p_1^{min(a_1,b_1)}\times p_2^{min(a_2,b_2)}\times ……p_1^{min(a_k,b_k)}\\ &lcm(a,b)=p_1^{max(a_1,b_1)}\times p_2^{max(a_2,b_2)}\times ……p_1^{max(a_k,b_k)}\\ &因此,gcd(a,b)\times lcm(a,b)=a\times b \end{align} \]
LL lcm(LL a, LL b) {
return a / gcd(a, b) * b;
}
例題:hdu1788
題意是求最小的正整數N,滿足除以每一個給定的數Mi,其結果均為Mi-a。作如下轉化:
\[ N\%M_i\equiv(M_i-a)\%M_i\Rightarrow (N+a)\%M_i=0(a<M_i<100) \]
即N是Mi的倍數,那麼題意也就是求所有M的最小公倍數,注意数字範圍即可。
同余
對於兩個不同的數a,b,如果有a % p = b % p(p>1),那我們就稱a和b模p同余,記作:
\[ a\equiv b(mod\space p) \]
同余的性質
\[ \begin{aligned} &1.自反性:a\equiv a(mod\space m)\\ &2.對稱性:a\equiv b(mod\space m),則b\equiv a(mod\space m)\\ &3.傳遞性:若a\equiv b(mod\space m),b\equiv c(mod\space m),則a\equiv c(mod\space m)\\ &4.線性運算:若a\equiv b(mod\space m),c\equiv d(mod\space m)\\ &\space\space\space則:a\pm c\equiv b\pm d(mod\space m);\space a\times c\equiv b\times d(mod\space m)\\ &5.冪運算:若a\equiv b(mod\space m),那麼a^n\equiv b^n(mod\space m)\\ &6.若ac\equiv bc(mod\space m),且c≠0,那麼a\equiv b(mod\space {m\over gcd(c,m)})\\ &\space\space\space特殊地,若gcd(c,m)=1,則a\equiv b(mod\space m) \end{aligned} \]
EXGCD
貝祖定理(Bezouts Identity):若設a,b是不全為0的整數,則存在整數x,y,使得ax+by=gcd(a,b),(a,b)代表最大公因數,則設a,b是不全為零的整數,則肯定存在整數x,y,使得ax+by=(a,b)。這個式子稱為貝祖等式。
EXtend GCD 即擴展歐幾里得算法,它可以用於解關於x,y的貝祖等式。
EXGCD的可行性
當a=0時,ax+by=gcd(a,b)=gcd(0,b)=b,這時可以解出y=1而x為任意。同理可得b=0時,解得x=1而y為任意。當a,b都不為0時:
\[ \begin{aligned} &設存在一組解為x_1,y_1,即ax_1+by_1=gcd(a,b)\\ &由歐幾里得定理:gcd(a,b)=gcd(b,a\%b)得:\\ &ax_1+by_1=gcd(a,b)=gcd(b,a\%b)\\ &那麼,bx_2+(a\%b)y_2=gcd(b,a\%b)=ax_1+by_1\\ &又a\%b=a-a/b\times b\\ &\begin {align} 那麼,&ax_1+by_1\\ =&bx_2+(a-a/b\times b)y_2\\ =&bx_2+ay_2-a/b*b*y_2\\ =&ay_2+b(x_2-a/b*y_2) \end{align}\\ &故可以得到:x_1=y_2,\space y_1=x_2-a/b*y_2 \end{aligned} \]
a或b為0即迭代求解的出口,迭代過程用代碼錶示如下:
typedef long long LL;
// ver 1
// 返回值為gcd(a,b),x和y即求出的特解
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(b, a % b, x, y);
//這裏通過迭代求出了一組解x,y,這組解對應上面推導過程中的x2和y2。
LL temp = x;
x = y;
y = temp - a / b * y;
return ans;
}
// ver 2
// 在熟悉了推導過程之後,我們可以利用引用的性質得到更為簡潔的ver2寫法
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(b, a % b, y, x);
y -= a / b * x;
return ans;
}
EXGCD的應用
求關於x,y的方程ax+by=c的一組解
這裏只需要對貝祖等式稍作變換即可。假設EXGCD求出的方程ax+by=gcd(a,b)的一組解為x1和y1。在方程兩邊同時乘上一個數m,使得c=m*gcd(a,b),這裏也就要求c是gcd(a,b)的倍數,即c%gcd(a,b)=0。這也是方程有解的條件。此時方程為:
\[ \begin{align} &ax_1\times m+by_1\times m=gcd(a,b)\times m=c \end{align} \]
對於ax+by=gcd(a,b),其參數解為:
\[ \begin{align} &x=x_1+{b\over gcd(a,b)}\times t,y=y_1-{a\over gcd(a,b)}\times t,t為參數\\ &這裏選用{b\over gcd(a,b)}和{a\over gcd(a,b)}是為了保證結果均為整數 \end{align} \]
那麼對比方程ax+by=c,可以得到其特解為:
\[ \begin{align} &x_0=x_1\times m=x_1\times {c\over gcd(a,b)}\\ &y_0=y_1\times m=y_1\times {c\over gcd(a,b)} \end{align} \]
那麼其通解為:
\[ \begin{align} &X=x_0+{b\over gcd(a,b)}\times t,Y=y_0-{a\over gcd(a,b)}\times t,t為參數\\ &這裏選用{b\over gcd(a,b)}和{a\over gcd(a,b)}作為係數是為了確保得到均為整數解。 \end{align} \]
對於任意一個確定的t,都有一組確定的解與之對應,只需要根據需要找出對應的解即可。
例如求滿足ax+by=c的X的最小非負整數解,那麼在得到X的通解之後只需調整t,調整過程用偽代碼錶示為:
S = b / gcd(a, b);
X = (x0 % S + S) % S;
// 由於x0可能是負數,故在模之後還要加上S在取一次模。
// 同理可得Y的最小非負整數解
T = a / gcd(a, b);
Y = (y0 % T + T) % T;
完整的exgcd求解方程ax+by=c的X和Y的最小非負整數解得代碼如下:
LL exgcd(LL a, LL b, LL c, LL& x, LL& y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(a, b, c, y, x);
y -= a / b * x;
LL S = b / ans, T = a / ans;
x = (x % S + S) % S;
y = (y % T + T) % T;
return ans;
}
求關於x的方程ax≡b(mod m)的最小非負整數解
ax≡b(mod m)即ax%m=b%m,即求ax+my=b%m(取模其實就相當於減去(加上)了若干個模數)。但只有a和m互質時有唯一解,否則無解。轉化為ax+my=b(假設b<m)后,套用exgcd即可。
求a關於模數p的乘法逆元
設x是a關於p的乘法逆元,那麼有ax≡1(mod p),這是第二個應用的特殊情況,顯然也可以通過類似方法求得。
逆元
Inverse Element,逆元,推廣了加法中的加法逆元和乘法中的倒數。直觀地說,它是一個可以取消另一給定元素運算的元素。a關於模p的逆元存在的條件是gcd(a,p)=1。
\[ \begin{align} &在模p意義下,設A的逆元為A^{-1},那麼有A\times A^{-1}\equiv 1(mod\space p) \end{align} \]
為什麼需要逆元呢?
\[ \begin{align} &在模意義下,{A\over B}\%p = {A\%p\over B\%p}\%p並不成立,\\ &那麼設B在模p意義下的逆元表示為B^{-1},根據逆元的定義,有{A\over B}=A\times B^{-1}(mod \space p)。 \end{align} \]
這裏,我們把除法轉化為乘法,就可以運用取模運算的性質:(a * b) % c = (a % c * b % c) % c,優化算法。
逆元的求法
EXGCD求法
給定模數m,求a的逆元相當於求解關於x的方程ax≡1(mod m)。這個方程可以轉化為ax-my=1 。用EXGCD可以求得一組x,y和gcd。檢查gcd是否為1(因為EXGCD是把問題轉化為求解方程ax-my=gcd(a,m),顯然只有gcd(a,m)=1時才能轉化),gcd不為1則說明逆元不存在。在能夠成功求解的情況下,把解調整x到0~m-1的範圍中即可。
費馬小定理
在模p為質數的情況下,有費馬小定理:
\[ a^{p-1}\equiv 1(mod\space p) \]
那麼:
\[ a\times a^{p-2}\equiv 1(mod\space p)\\ 故a的逆元a^{-1}=a^{p-2} \]
然後用快速冪求出逆元即可。
拓展:快速冪
// a是底數,b是指數
LL pow(LL a, LL b, LL mod) {
LL ans = 1;
while(b) {
if(b&1) ans = ans * a % mod;
b >>= 1, a = a * a % mod;
}
return ans;
}
拓展:慢速乘,防止快速冪過程中乘法運算爆long long,只需把快速冪中的乘法替換為mul()即可,一般情況用不上。
// a為因數1,b為因數2
LL mul(LL a, LL b, LL mod) {
LL ans = 0;
while(b) {
if(b&1) ans = (ans + a) % mod;
b >>= 1, a = (a + a) % mod;
}
return ans;
}
歐拉定理
費馬小定理其實是歐拉定理的特殊情況。在模數p不為質數時,有:
\[ a^{\phi(p)}\equiv 1(其中\phi()是歐拉函數) \]
同理可以得到:
\[ a^{-1}\equiv a^{\phi(p)-1} \]
線性逆元表
有時我們需要快速得出1~p-1的所有逆元,這個時候我們需要一種O(n)的方法,這就是線性逆元表。
首先,1關於任意模p的逆元的逆元都是1。設p = k * i + r(r < i ,1 < i < p),那麼有:
\[ k\times i + r \equiv0(mod\space p) \]
利用逆元性質,兩邊同時乘上i的逆元以及r的逆元,得到:
\[ k\times r^{-1}+i^{-1}\equiv 0(mod\space p) \]
移項得:
\[ i^{-1}\equiv-k \times r^{-1}\equiv -⌊{p\over i}⌋\times (p\%i)^{-1}(mod\space p) \]
顯然p%i是小於i的,那麼我們就可以利用之前的結果在O(1)時間算出單獨的逆元。
用代碼錶示即:
// inv[i]對應i的逆元
LL inv[maxn];
void CalInv() {
// 0沒有逆元,故不初始化
inv[1] = 1;
for (int i = 2; i < maxn; i++)
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
利用
\[ (a\times c)\% mod = (a\%mod\times c\%mod)\%mod;\\ 由a\times a^{-1}\equiv1(mod \space p)得,\\ inv_{fac(i)}\equiv inv_{fac(i)}\times (i+1)^{-1}\times (i+1)\equiv inv_{fac(i+1)}\times (i+1)(mod\space p) \]
我們還可以在線性時間求出1~min(n,p)的階乘的逆元,代碼如下:
LL fac[maxn], inv[maxn];
void CalFacInv() {
// 0的階乘是1
fac[0] = fac[1] = 1;
for (int i = 2; i < maxn; i++)
fac[i] = fac[i - 1] * i % mod;
inv[maxn - 1] = QPow(fac[maxn - 1], mod - 2);
// 注意邊界是≥
for (int i = maxn - 2; i >= 0; i--)
inv[i] = inv[i + 1] * (i + 1) % mod;
}
例題:hdu1576
題意即把除法轉化為逆元乘法。3種方法的代碼:
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
LL n, B;
const LL mod = 9973;
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
LL ans = exgcd(b, a % b, y, x);
y -= a/b*x;
return ans;
}
LL QPow(LL x, LL n)
{
LL res = 1;
while(n)
{
if(n & 1) res = res * x % mod;
x = x * x % mod;
n >>= 1;
}
return res;
}
LL inv[10000];
void CalInv() {
inv[1] = 1;
for (LL i = 2; i < 10000; i++) {
inv[i] = (mod * mod - mod / i * inv[mod % i]) % mod;
}
}
int main() {
CalInv();
int t;
cin >> t;
while (t--) {
cin >> n >> B;
// exgcd
// LL x, y;
// exgcd(B, mod, x, y);
// cout << (x + mod) * n % mod << endl;
// QPow
// cout << QPow(B, mod - 2) * n % mod << endl;
// 逆元打表
cout << inv[B % mod] * n % mod << endl;
}
return 0;
}
本站聲明:網站內容來源於博客園,如有侵權,請聯繫我們,我們將及時處理
【其他文章推薦】
※USB CONNECTOR掌控什麼技術要點? 帶您認識其相關發展及效能
※評比前十大台北網頁設計、台北網站設計公司知名案例作品心得分享
※智慧手機時代的來臨,RWD網頁設計已成為網頁設計推薦首選
※台灣海運大陸貨務運送流程
※兩岸物流進出口一站式服務