初等數論——素數，逆元，EXGCD有關

初等數論

素數定義

設整數 \(p\ne 0,\pm 1\) 。如果 \(p\) 除了平凡約數以外沒有其他約數，那麼稱 \(p\) 為素數（不可約數）。

若整數 \(a\ne 0,\pm 1\) 且 \(a\) 不是素數，則稱 \(a\) 為合數。

——————OI Wiki

素數的判定（素性測試）

如何判斷一個數 \(x\) 是不是質數？

很顯然我們可以暴力的列舉 \(1\) 到 \(\sqrt{x}\) 來看是否整除來判斷，但複雜度太高了。

#include<bits/stdc++.h>
#define int long long 
using namespace std;
int t,n;
inline int pd(int x)
{
    int len=sqrt(x);
    for(int i=2;i<=len;i++)
        if(x%i==0)return 0;
    return 1;
}
signed main()
{
    cin>>t;
    while(t--)
    {
        cin>>n;
        if(pd(n))cout<<"Y"<<endl;
        else cout<<"N"<<endl;
    }
    return 0;
}

Miller-Rabin素性測試

這個演算法就一個用法：判斷一個數是不是素數。

這個演算法的依據是：

費馬小定理

對於任意的質數 \(p\),任意整數 \(a\),均滿足:\(a^{p}\equiv a\pmod{p}\)，如果 \(a\) 不是 \(p\) 的倍數，這個也可以寫成 \(a^{p-1} \equiv 1 \pmod{p}\)

證明如下:
這個式子可以用尤拉定理來證明:

首先,把式子稍微變換一下得到:

\(a^{p-1}\times a\equiv a\pmod{p}\),因為 \(a\equiv a\pmod{p}\) 恆成立,所以 \(a^{p-1}\bmod p=1\) 時費馬小定理才成立,又因為 \(p\) 是質數,所以 \(\varphi(n)=n-1\),所以根據尤拉定理:若 \(a\),\(p\) 互質則 \(a^{p-1}\bmod p=1\) 成立。那麼對於 \(a\),\(p\) 不互質，因為 \(p\) 是質數，所以，\(a\) 一定是倍數 \(a^{p}\equiv a\equiv 0 \pmod{p}\)。綜上所述，費馬小定理得證，其實算是一個尤拉定理的特例。

那麼是不是當一個數 \(p\) 滿足任意 \(a\) 使得 \(a^{p-1}\equiv \pmod{p}\) 成立的時候它就是素數呢？

答案是 false

所以我們需要下面的定理來減小錯誤的機率。

二次探測定理

若 \(p\) 為素數，\(a^{2}\equiv \pmod{p}\)，那麼 \(a\equiv \pm 1\pmod{p}\)

證明：

\(a^{2}\equiv 1\pmod{p}\)

\(a^{2}-1\equiv 0\pmod{p}\)

\((a-1)\times (a+1)\equiv 0\pmod{p}\)

\((a-1)\equiv 0\pmod{p}\) 或 \((a+1)\equiv 0\pmod{p}\)

\(a\equiv \pm 1\pmod{p}\)

演算法過程

我們把 \(p-1\) 分解成 \(2^{k}\times t\) 的形式

能這麼分的原因是 \(p\) 如果是素數，那他只要不是 \(2\)，那麼 \(p-1\) 一定是偶數，也就是都至少能分離出一個二出來

當 \(p\) 是素數，有 \(a^{2^{k}\times t}\equiv 1\pmod{p}\)

然後隨機選擇一個數 \(a\)，計算出 \(a^{t}\pmod{p}\)

然後如果要是當前的值就是 \(1\)，那麼說明後面再平方也還是 \(1\)（二次探測定理），可以直接返回判斷是素數；否則就開始跟據計算出來的 \(k\) 來進行 \(k\) 次的平方操作，如果要是遇到值為 \(p-1\) 的值，下一次平方值必定為 \(1\)，所以也可以直接返回，這樣就實現了利用費馬小定理的同時利用二次探測定理來判斷素數了

百度告訴我們若 \(p\) 透過一次測試，則 \(p\) 不是素數的機率為 \(25\%\)

那麼我們用 \(2,3,5,7,11,13,17,19,23,27...\) 來進行多次判斷，正確性就會大大提高。

當然也不會是 \(100\%\)，但是有以下的結論。

我們需要注意的是，因為是次方運算有點猛，所以推薦使用龜速乘來防止溢位。

埃氏篩

從小到大列舉範圍，如果當前數沒有被標記就說明是質數，直接加進去，然後列舉在範圍內的當前數的倍數標記，然後一直重複直到結束，算是比較優秀的篩了（至少相對暴力）據說複雜度 \(O(n\log\log n)\)

圖片來自 b 站董曉演算法。

線性篩

重中之重，必須掌握——————我說的

當然我覺得這個還是更快一些的。

我們首先要知道的就是上面的篩慢在什麼地方？

因為有很多數被重複篩了好幾次，那麼我們有什麼辦法能不能讓一個數只被篩一次，讓複雜度成為真正的線性呢？

答案是 true

我們和上面一樣，只不過在篩的時候，我們不是列舉當前質數的多少倍了，而是每列舉到一個數 \(i\)，然後就開始列舉當前篩出來的質數的 \(i\) 倍，而噹噹前 \(i\) 對當前列舉到的質數取模是 \(0\) 的時候，就直接退出，這樣可以做到每一個合數都被自己的最小質因數篩掉，因為當前 \(i%p=0\) 說明 \(i\) 已經被 \(p\)，也就是 \(i\) 的最小質因數給篩掉了，那麼後面在出現的 \(i\) 乘以各個素數的數一定會被 \(p\) 給篩掉，現在就不用篩了，直接退出。

code:

#include<bits/stdc++.h>
using namespace std;
int n,m,ans[100010000],vis[100010000],cnt;
int main()
{
	std::ios::sync_with_stdio(0);
	cin>>n>>m;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])ans[++cnt]=i;
		for(int j=1;j<=cnt&&ans[j]*i<=n;j++)
		{
			vis[ans[j]*i]=1;
			if(i%ans[j]==0)break;
		}
	}
	for(int i=1;i<=m;i++)
	{
		int x;
		cin>>x;
		cout<<ans[x]<<endl;
	}
	return 0;
}

EXGCD

歐幾里得應該在上一篇關於數論的文章裡講到了，而且也不難，所以在此不再贅述。

基本演算法：對於不完全為 \(0\) 的非負整數 \(a\),\(b\),\(\gcd(a,b)\) 表示 \(a\),\(b\) 的最大公約數，必然存在整數對 \(x\),\(y\)，使得 \(\gcd(a,b)=a\times x+b\times y\).

證明：

設 \(a>b\).

顯然當 \(b=0,\gcd(a,b)=a\),此時 \(x=1,y=0\).

當 \(a\times b \ne 0\) 時；

設 \(a\times x_{1}+b\times y_{1}=\gcd(a,b)\);

\(b\times x_{2}+(a\bmod b)\times y_{2}=\gcd(b,a\bmod b)\)

根據樸素的歐幾里得原理有 \(\gcd(a,b)=gcd(b,a\bmod b)\);

則：\(a\times x_{1}+b\times y_{1}=b\times x_{2}+(a\bmod b)\times y_{2}\)

即：\(a\times x_{1}+b\times y_{1}=b\times x_{2}+(a-\left \lfloor\frac{a}{b}\right \rfloor\times b)\times y_{2}=a\times y_{2}+b\times x_{2}-\left \lfloor\frac{a}{b}\right \rfloor\times b\times y_{2}\)

等量代換得：\(x_{1}=y_{2};y_{1}=x_{2}-\left \lfloor\frac{a}{b}\right \rfloor\times y_{2}\);

這樣我們就可以一直遞迴求 \(x1,y1\) 了。

因為總有一個時候 \(b=0\)，所以遞迴可以結束。

#include<bits/stdc++.h>
#define int long long 
using namespace std;
int x,y,n,m;
inline void exgcd(int a,int b)
{
	if(b==0){x=1,y=0;return ;}
	exgcd(b,a%b);
	int t=x;
	x=y;
	y=t-a/b*y;
}
signed main()
{
	cin>>n>>m;
	exgcd(n,m);
	x=(x%m+m)%m;
	cout<<x<<endl;
	return 0;
}

或者這種寫法：

#include<bits/stdc++.h>
#define int long long 
using namespace std;
int x,y,n,m;
inline void exgcd(int a,int b,int &x,int &y)
{
	if(b==0){x=1,y=0;return ;}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
signed main()
{
	cin>>n>>m;
	exgcd(n,m,x,y);
	x=(x%m+m)%m;
	cout<<x<<endl;
	return 0;
}

逆元

主要還是模運算下乘法運算的逆元，下文中逆元都會帶有下標 \(inv\)。

如果一個線性同餘方程 \(ax\equiv 1\pmod{b}\)，則稱 \(x\) 為 \(a\bmod b\) 的逆元，記作 \(a^{-1}\)。
——————OI Wiki

EXGCD求逆元

\(a \equiv b \pmod{m}\)

老朋友了，表示 \(a\bmod m=b\mod m\)。

然後看一下這個：\(\frac{a}{b} \equiv x\pmod{p}\).

模意義下的除法在大多數時候都不適用。

當題目中要求對答案取模，但我們又不得不使用一般的除法的時候，就需要用逆元的乘法來代替除法。

根據逆元的定義得：

\((a\times b_{inv}) \equiv x \pmod{p}\);

左右兩邊同乘 \(b\),移項得：

\(\left\{\begin{matrix} (b\times a\times b_{inv})\equiv x\times b \pmod{p}\\ a\equiv b\times x \pmod{p} \end{matrix}\right.\)

上減下得：

\(a\times(b\times b_{inv}-1)\equiv 0 \pmod{p}\)

移項得：

\(b\times b_{inv}-1\equiv 0 \pmod{p}\);

兩邊同時加一得：

$ b\times b_{inv}\equiv 1\pmod{p}$；

設一個正整數 \(k\) ；

兩邊同時加 \(k\times p\) 得；

$ b\times b_{inv}+k\times p\equiv 1+k\times p\pmod{p}$；

因為 \((k\times p )\bmod p=0\);

所以：$ b\times b_{inv}+k\times p\equiv 1\pmod{p}$；

根據我們的擴充套件歐幾里得可以得知：當知道 \(a,b\) 的值的時候是可以求出 \(a\times x+b\times y=gcd(a,b)\) 中的 \(x,y\) 的值。當然前提是 \(b,p\) 互質的時候上面的柿子就可以得出：

\(b\times b_{inv}+k\times p= 1\)；

知道 \(b,p\) 就可以愉快的求 \(b_{inv}\) 啦。

好的到了這一步已經可以有兩種解決方式了，一種是用擴充套件歐幾里得來解逆元，另一種是用費馬小定理，先來好理解的。

費馬小定理求逆元

首先回到我們前面推的柿子；

$ b\times b_{inv}\equiv 1\pmod{p}$；

好的現在我們可以知道當 \(p\) 為質數時：

\(a^{p}\equiv a\pmod{p}\);

稍微變換一下得：

\(a^{p-1}\times a\equiv a\pmod{p}\)

因為 \(a\) 是正整數，兩邊同除 \(a\) 得：

\(a^{p-1}\equiv 1 \pmod{p}\);

根據我們上面得出得柿子一結合就可以得出：

\(\left.\begin{matrix} b\times b_{inv}\equiv 1\pmod{p} \\ b\times b^{p-2}\equiv 1\pmod{p} \end{matrix}\right\}\Rightarrow b_{inv}=b^{p-2}\)

所以我們就可以用：

\(a\times b^{p-2} \bmod p\) 來求 $ \frac{a}{b}\bmod p$ 啦。

線性求逆元

前面的都太慢了，我們要追求效率。

首先我們有一個 \(1_{inv}\equiv 1\pmod{p}\);

我們設 \(p=k\times i+r\)；

然後放到 \(\pmod{p}\) 意義下就可以得到：

\(k\times i+r\equiv 0\pmod{p}\);

乘上 \(i_{inv},r_{inv}\) 可以得到;

\(k\times r_{inv}+i_{inv}\equiv 0\pmod{p}\)

移項得：

\(i_{inv}\equiv -k\times r_{inv}\pmod{p}\)

把 \(k,r_{inv}\) 用 \(p\) 表示出來就好了。

\(i_{inv}\equiv -\left \lfloor \frac{p}{i} \right \rfloor \times (p\bmod i)_{inv} \pmod{p}\);

貼一下程式碼：

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

積性函式

數論函式：在數論上，對於定義域為正整數，值域為複數的函式，我們稱之為數論函式。

非完全積性函式：對於數論函式 \(f\)，若滿足 \(gcd(a,b)=1\) 時，有 \(f(ab)=f(a)\times f(b)\)，則稱函式 \(f\) 為積性函式.

完全積性函式：對於數論函式 \(f\)，即使不滿足 \(gcd(a,b)=1\)，仍有 \(f(ab)=f(a)\times f(b)\)，則稱函式 \(f\) 為完全積性函式.

簡單積性函式

\(\varphi (n)\)————尤拉函式，計算從一到 \(n\) 與 \(n\) 互質的數的數目。

\(\mu(n)\)————莫比烏斯函式，關於非平方數的質因子數目.

\(\mu(n)=\left\{\begin{matrix} 1 \ \ (n=1) \\ 0\ \ \ \ (n\text{含有1除外的平方因子}) \\ (-1)^{k}\ \ \ \ (k\text{為}n\text{的本質不同質因子個數）} \end{matrix}\right.\)

\(\gcd(n,k)\)————最大公約數函式，當 \(k\) 固定的情況時是積性函式。

\(d(n)\)————\(n\) 的正因子數目。

\(\sigma(n)\)————\(n\) 的正因子之和。

\(\sigma k(n)\)————因子函式，\(n\) 的所有正因子的 \(k\) 次冪之和，當中 \(k\) 可以為任何複數。

\(1(n)\)————不變的函式，定義為 \(1(n)=1\) （完全積性函式）。

\(id(n)\)————單位函式，定義為 \(id(n)=n\) （完全積性函式）。

\(idk(n)\)————冪函式，對於任何複數、實數 \(k\)，定義為 \(idk(n)=n^k\)（完全積性函式）。

\(\xi(n)\)————定義為：若 \(n=1\),\(\xi(n)=1\);若 \(n>1,\xi(n)=0\)。別稱為“對於狄利克雷卷積的乘法單位”（完全積性函式）。

\(\lambda(n)\)————劉維爾函式，關於能整除 \(n\) 的質因子的數目.

\(\gamma(n)\)————定義為 \(\gamma(n)=(-1)^{\omega(n)}\),在此加性函式 \(\omega(n)\) 是不同能整除 \(n\) 的質數的數目。

參考資料：OI WIKI，自為風月馬前卒的部落格，https://www.cnblogs.com/zylAK/p/9569668.html ，b站董曉演算法。