1. 程式人生 > >關於數論乘法逆元及相關知識點

關於數論乘法逆元及相關知識點

在求解a/b%m時,可以轉化為(a%(b*m))/b,轉化過程如下

令k = (a/b)/m(向下取整), x = (a/b)%m;

a/b = k*m + x (x < m);

a = k*b*m + b*x;

a%(b*m) = b*x;

a%(b*m)/b = x;

得證: a/b%m = a%(b*m)/b;(公式適用於很多情況:m不必是素數,b和m也不必互質)

上面的公式適用於b較小,a需要線上求且較大的時候;

比如:求a/b的值,且另a = x^k, 求a的過程隨時會爆精度所以需要對a進行取模,但是求的過程中如果是直接對a進行%mod,之後再進行a/b,會出現錯誤的。正確的做法是根據上述公式,在求a過程中對(b*m)取模,之後再除以一次b即可。[ a/b%m = (a%(b*m))/b%m ]

所有的除法取模問題都可以用這種方法,但是當b很大的時候,則會出現爆精度問題,所以引出乘法逆元,將除法取模轉換為乘法取模。

b存在乘法逆元的充要條件是b與模數m互質。設令c為b的逆元,即b*c ≡ 1(mod m);

解釋:滿足同餘方程,b*c和1對m的模數相同,即b*c對m取餘為1(m > 1);

那麼 a/b = (a/b)*1 = (a/b)*(b*c)(mod m) = a*c(mod m);

即,除以一個數對m取模等於乘以這個數的逆元對m取模;

三種求逆元的方法&:

1.逆元求解利用擴歐。
2.當m為質數的時候直接使用費馬小定理,m非質數時使用尤拉函式。
3.當m為質數的時候,使用神奇的線性方法。

擴充套件&:

1.利用乘法逆元求解組合數及組合數的其它求法。

2.線性時間求階乘的逆元。

擴充套件歐幾里得演算法.

首先了解一下擴充套件歐幾里得演算法,

已知 a 和 b 的最大公約數是 gcd,那麼,我們一定能夠找到這樣的 x 和 y ,使得: a*x + b*y = gcd 這是一個不定方程(其實是一種丟番圖方程),有多解是一定的,但是隻要我們找到一組特殊的解 x0 和 y0 那麼,我們就可以用 x0 和 y0 表示出整個不定方程的通解:
    x = x0 + (b/gcd)*t; ((b/gcd)是b的因子,能表示的數只能多而不會少,下面同理)
    y = y0 – (a/gcd)*t;

現在,我們知道了一定存在 x 和 y 使得 : a*x + b*y = gcd , 那麼,怎麼求出這個特解 x 和 y 呢?只需要在歐幾里德演算法的基礎上加點改動就行了。
我們觀察到:歐幾里德演算法停止的狀態是: a = gcd , b = 0 ,那麼,這是否能給我們求解 x, y 提供一種思路呢?這時候,只要 a = gcd 的係數是 1 ,只要 b 的係數是 0 或者其他值(無所謂是多少,反正任何數乘以 0 都等於 0 但是a 的係數一定要是 1),這時,我們就會有: a*1 + b*0 = gcd;
當然這是最終狀態,但是我們是否可以從最終狀態反推到最初的狀態呢?
假設當前我們要處理的是求出 a 和 b的最大公約數,並求出 x 和 y 使得 a*x + b*y= gcd ,而我們已經求出了下一個狀態:b 和 a%b 的最大公約數,並且求出了一組x1 和y1 使得: b*x1 + (a%b)*y1 = gcd , 那麼這兩個相鄰的狀態之間是否存在一種關係呢?
我們知道: a%b = a - (a/b)*b(這裡的 “/” 指的是整除,例如 5/2=2 , 1/3=0),那麼,我們可以進一步得到:
        gcd = b*x1 + (a-(a/b)*b)*y1;
            = b*x1 + a*y1 – (a/b)*b*y1;
            = a*y1 + b*(x1 – a/b*y1);
對比之前我們的狀態:求一組 x 和 y 使得:a*x + b*y = gcd,是否發現了什麼?
這裡:
        x = y1;
        y = x1 – a/b*y1;

上程式碼:

遞迴版e_gcd();


還有非遞迴版的;


前提條件:b與m互質,存在唯一解。

擴充套件歐幾里德演算法,我們既可以求出最大公約數,還可以順帶求解出使得: a*x + b*y = gcd 的通解 x 和 y
擴充套件歐幾里德有什麼用處呢?
求解形如 a*x +b*y = c 的通解,但是一般沒有誰會無聊到讓你寫出一串通解出來,都是讓你在通解中選出一些特殊的解,乘法逆元便是。

a*x ≡ 1(mod m);

這裡稱 x 是 a 關於 m 的乘法逆元, 這怎麼求?可以等價於這樣的表示式: a*x + m*y = 1;
看出什麼來了嗎?沒錯,當gcd(a , m) != 1 的時候是沒有解的, 這也正是 a*x + b*y = c 有解的充要條件: c % gcd(a , b) == 0;
接著乘法逆元講,一般,我們能夠找到無陣列解滿足條件,但是一般是讓你求解出最小的那組解,怎麼做?我們求解出來了一個特殊的解 x0 那麼,我們用 x0 % m其實就得到了最小的解了。為什麼?
可以這樣思考:
x 的通解不是 x0 + m*t 嗎?
那麼,也就是說, a是和 a 關於 m 的逆元是關於 m 同餘的,那麼根據最小整數原理,一定存在一個最小的正整數,它是 a 關於m 的逆元,而最小的肯定是在(0 , m)之間的,而且只有一個,這就好解釋了。
有時候我們得到的特解 x0 是一個負數,還有的時候我們的 m 也是一個負數這怎麼辦?
當 m 是負數的時候,我們取 m 的絕對值就行了,當 x0 是負數的時候,他模上 m 的結果仍然是負數(在計算機計算的結果上是這樣的,雖然定義的時候不是這樣的),這時候,我們仍然讓 x0 對abs(m) 取模,然後結果再加上 abs(m) 就行了,於是,我們不難寫出下面的程式碼求解一個數 a 對於另一個數 m 的乘法逆元。


當然略微有點長了,比賽時用下面的就足以解決很多。

int ex_gcd(int a, int b, int &d, int &x, int &y){
	if(b == 0){
		x = 1; y = 0; d = a;
		return;
	}
	ex_gcd(b, a%b, d, y, x);	//不斷交錯x和y
	y -= a/b*x;	//當遞迴完成return時,y值return給了x, x值return給了y
}
int clac(int a, int m){
	int x, y, d;
	ex_gcd(a, m, d, x, y);
	return (m + x%m)%m;
}

通過費馬小定理或者尤拉函式去求.

1.費馬小定理求逆元:

在m是素數的情況下,對任意整數x都有x^m ≡ x(mod m)。 
如果x為整數且無法被m整除,則有x^(m−1) ≡ 1(mod m)。 
所以可以在m為素數的情況下求出一個數的逆元,x * x^(m−2) ≡ 1(mod m),x^(m−2)即為逆元。

推導過程如下:

然後就可通過快速冪求得逆元。

quickM(int base, int b);

2.尤拉函式求逆元:
令ϕ(m)表示小於等於m且與m互素的正整數的個數(特例(ϕ(1) = 1))。
如果x和m互質,則有x^ϕ(m) ≡ 1(mod m),即x * x^(ϕ(m)−1) ≡ 1(mod m),x^(ϕ(m)−1)即為x的逆元。 
在m為質數的情況下,ϕ(m) = m−1,即為費馬小定理。

先介紹一下尤拉函式的一些知識(能學一點是一點...)

φ(x) = x*(1-1/p(1))(1-1/p(2))(1-1/p(3))(1-1/p(4))…..(1-1/p(n)) 其中p(1),p(2)…p(n)為x的所有質因數;x是正整數; 注意:每種質因數只有一個。φ(1)=1(唯一和1互質的數,且小於等於1)。

尤拉函式的性質:
(1)p^k型尤拉函式:
若N是質數p(即N=p),φ(N) = φ(p) = p-p^(1-1) = p-1。
若N是質數p的k次冪(即N=p^k),φ(N)=p^k-p^(k-1) = (p-1)*p^(k-1)。
(2)m*n型尤拉函式:
設n為正整數,以φ(n)表示不超過n且與n互素的正整數的個數,稱為n的尤拉函式值。若m,n互質,φ(m*n)=(m-1)*(n-1)=φ(m)*φ(n)。
(3)特殊性質:
若n為奇數時,φ(2n)=φ(n)。
(4)尤拉定理:
對於任何兩個互質的正整數x,m(m>2)有:x^φ(m) ≡ 1(mod m)。
(5)擴充套件(費馬小定理)
當m=p且x與素數p互質(即:gcd(x,p)=1)則上式有: x^(p-1) ≡ 1(mod p)。
(6)尤拉函式的延伸:
小於或等於n的數中,與n互質的數的總和為:φ(n) * n/2 (n>1)。

證明就不給出了,下面是模板程式碼:

求單個數的尤拉函式值:

int Euler(int n){
	int ans = n;
	//直接掃一般會超時,所以這裡利用任何一個合數都至少
	//有一個不大於根號n的素因子,所以只需遍歷到根號n即可 
	for(int i = 2; i*i <= n; ++i){
		if(n%i == 0){
			ans = ans/i*(i-1);	//先除後乘防止中間資料溢位
			while(n%i == 0) n /= i; 
		}
	}
	if(n > 1) ans = ans/n*(n-1);
	return ans;
}
多次用到尤拉函式值,可通過篩選預先打表(時間複雜度O(n·logn)):
void getEuler(){
	int euler[maxn];
	for(int i = 1; i <= maxn; ++i)
		euler[i] = i;
	for(int i = 2; i <= maxn; i+=2)
        euler[i] /= 2; 	//其實相當於進行了一次euler[i]/2*(2-1);
	for(int i = 3; i <= maxn; ++i)
    {
        if(euler[i] == i) //未被篩到, 即是素數, 則用此素數來篩
        {
            for(int j = i; j <= maxn; j+=i)
            {
                euler[j] = euler[j]/i*(i-1);
            }
        }
    }
} 

Last,

求出模數m的尤拉函式值後,進行一次快速冪 quickM(base,φ(n)-1)求得逆元。

線性時間內求範圍內所有整數的逆元.

當m不是很大的時候(陣列能夠存起來的時候)可以線性時間求逆元。規定m為質數,
首先 1^(−1) ≡ 1(mod m)
然後我們設 m = k*x + r, r < x, 1 < x < p;
再將這個式子放到mod m意義下就會得到
k*x + r ≡ 0 (mod m)
兩邊同時乘上 x^(−1) * r^(−1)就會得到
k * r^(−1) + x^(−1) ≡ 0(mod m)
x^(−1) ≡ −k * r^(−1)(mod m)
x^(−1) ≡ -(m/x) * (m%x)^(−1)(mod m)
從頭開始掃一遍即可,時間複雜度O(n);

int inv[maxn];
inv[1] = 1;
for(int i = 2; i < maxn; ++i)
    inv[i] = (-p/i + p) * inv[p%i] % p;

擴充套件:

利用乘法逆元求解組合數及組合數的其它求法.

首先宣告C(n, m)是指n箇中取m個,網上都沒統一寫法,真是的...

求組合數常用的公式為C(n, m) = n!/((n-m)! * m!);

當n, m都很小的時候可以利用楊輝三角直接求。
或者通過下面推導公式遞迴求。
C(n+1, m) = C(n, m) + C(n, m-1);
C(n, m) = C(n-1, m) + C(n-1, m-1);

簡單的遞迴程式碼:

#include <stdio.h>
int calc(int n, int m){
	if(n < m || n == 0) return 0;
	if(n == m || m == 0) return 1;
	return calc(n-1, m) + calc(n-1, m-1);
}
int main(){
	int n, m;
	while(scanf("%d %d", &n, &m)){
		printf("%d\n", calc(n, m));
	}
	return 0;
}

再就是最常用的利用乘法逆元進行求解。

//可預先打出指定範圍內階乘的表
void init(){
	fact[0] = 1;
	for(int i = 1; i <= maxn; ++i)
	fact[i] = fact[i-1]*i % mod;
}

1. (n!/(m!*(n-m)!)) % mod = x %mod ,先對算出n!、m!、(n-m)!對mod取模的餘數,就轉換為a/b = x%mod;因為m為素數,所以等價於b*x +mod*y = gcd(b,mod); 然後用擴充套件的歐幾里得定理算出 b*x0 +mod*y0 = 1的特解x0,x再進而確定 =(mod + x0%mod) %mod; 則a/b = a*x%mod;

LL e_gcd(LL a, LL b, LL &x, LL &y){
	if(b == 0){
		x = 1; y = 0;
		return a;
	}
	LL ans = e_gcd(b, a%b, y, x);
	y -= x*(a/b);
	return ans;
}
LL clac(LL a, LL m){
	LL x, y;
	e_gcd(a, m, x, y);
	return (m + x%m)%m;
}
LL C(int n, int m){
	return fact[n]*clac(fact[m]*fact[n-m]%mod, mod)%mod;
}

2.如果mod是素數 則b的逆元其實就是b^(mod-2)即 (m!*(n-m)!)的逆元為 (m!*(n-m)!)^(mod-2);

LL quickM(LL a, LL b){
	LL base = a%mod, ans = 1;
	while(b){
		if(b&1) ans = ans*base%mod;
		base = base*base%mod;
		b >>= 1; 
	}
	return ans;
}
LL C(int n, int m){
	return fact[n]*quickM(fact[m]*fact[n-m]%mod, mod-2)%mod;
}

3.當n和m比較大而mod為素數且比較小(10^5左右)的時候,可以用Lucas定理計算.

Lucas定理:

A、B是非負整數,模數mod = p(p是質數)。A、B寫成p進位制:
A = a[n]a[n-1]...a[0];
B = b[n]b[n-1]...b[0];
則組合數C(A,B)與C(a[n], b[n])*C(a[n-1], b[n-1])*...*C(a[0], b[0])(mod p)同餘
即:Lucas(n, m, p) = c(n%p, m%p)*Lucas(n/p, m/p, p);

模板程式碼:

LL quickM(LL base, LL b, LL p) {
	LL ans = 1;
	while(b) {
	    if(b&1) ans = (ans * base) % p;
	    base = (base*base) % p;
	    b >>= 1;
	}
	return ans;
}

//n,m過大不能打表只能在線求階乘
LL Comb(LL a, LL b, LL p) {
	if(a < b) return 0;
	if(a == b) return 1;
	if(b > a - b) b = a - b;
	
	LL ans = 1, ca = 1, cb = 1;
	for(LL i = 0; i < b; ++i) {
	    ca = (ca * (a-i)) % p;
	    cb = (cb * (b-i)) % p;
	}
	ans = (ca*quickM(cb, p-2, p)) % p;
	return ans;
}

LL Lucas(int n, int m, int p) {
	LL ans = 1;
	while(n && m && ans) {
		ans = (ans*Comb(n%p, m%p, p)) % p;
		n /= p;
		m /= p;
	}
	return ans;
}

4.打表求階乘數逆元的新方法.

打一個1~n的階乘的逆元的表,假如n!的逆元叫做f[n],可以先用費馬小定理、擴充套件歐幾里得等 
求出f[n],再用遞推公式求出前面的項。
我們記數字 x 的逆元為 f(x) (mod m)。
  因為 n! = (n-1)! * n;
  所以 f(n!) = f((n-1)! * n) = f((n-1)!) * f(n);
  所以 f((n-1)!) = f(n!) * f(f(n)) = f(n!) * n;  (數的逆元的逆元就是它自身)
這樣子我們就可以用後項推出前面的項了。

LL quickM(LL base, LL b)
{  
    LL ans = 1;  
    while(b)
	{  
        if(b&1) ans = (ans * base) % mod;  
        base = (base*base) % mod;  
        b >>= 1;
    }  
    return ans; 
}
void init()
{
    fact[0] = 1;
    for(int i = 1; i <= maxn; ++i)
    fact[i] = fact[i-1]*i%mod;
    fiv[maxn] = quickM(fact[maxn], mod-2);
    for(int i = maxn-1; i >= 0; --i)
    {
        fiv[i] = fiv[i+1]*(i+1);
        fiv[i] %= mod;
    }
}
LL C(int n, int m)
{
	if(m > n) return 0;
	return fact[n]*fiv[n-m]%mod*fiv[m]%mod;
}

5.發現一個近乎完美的init(): 階乘、整數逆元、階乘逆元,一次O(n)即可。

void init()
{
	fact[0] = fact[1] = 1;
	fiv[0] = fiv[1] = 1; 
	inv[1] = 1;
	for(int i = 2; i <= maxn; ++i)  
	{
		//遞推儲存fact階乘,遞推求inv和fiv各個逆元 
	    fact[i] = fact[i-1]*i%mod;
	    inv[i] = (mod-mod/i)*inv[mod%i]%mod;
	    fiv[i] = inv[i]*fiv[i-1]%mod;
	}
}

Over~~~

參考博文: