數論專題

數論複習（一）
內容概要：素數篩法，尤拉函式，逆元，擴充套件中國剩餘定理，取模，拓展歐幾里得（質因數分解與輾轉相除），裴蜀定理，線性同餘方程，莫特蘭-切比雪夫定理，卡特蘭數
數論複習（二）
內容概要：逆元，素數篩法，素數判斷，質因數分解
背誦程式碼
內容概要：快速冪，擴充套件gcd，過載運算子【有關數論】

[數論複習（三）]
內容概要：組合數
數論複習（四）
內容概要：概率，期望

0.gcd

求最大公因數

法1：更相減損術

//下面這個gcd函式在正int型內完全通用，返回a,b的最大公因數。
//但是當a,b之間差距較大時（如100000倍）會導致錯誤（棧過深）
int gcd(int a,int b){
    if(a==b)return a;
    else if(a>b)a-=b;
    else b-=a;
    return gcd(a,b);
}
int main(){
    int a,b;
    cin>>a>>b;
    cout<<gcd(a,b)<<endl;
    return 0;
}
 

法2：輾轉相除法（歐幾里得演算法）

int gcd(int a,int b)
{
	if(b==0) return a;
	else return gcd(b,a%b);
}

直接在法1改進，效率倍增

法3：一個接近定理的東西：質因數分解

對\(a,b\)進行素因子分解：
\(a=p_1^{r_1}p_2^{r_2}...p_n^{r_n},b=p_1^{s_1}p_2^{s_2}...p_n^{s_n}\)

注意：r1，s2等代表質因子需要乘的次數

則有

\[\gcd (a,b)=p_1^{\min (r_1,s_1)}p_2^{\min (r_2,s_2)}...p_n^{\min (r_n,s_n)} \]

\[\text{lcm} (a,b)=p_1^{\max (r_1,s_1)}p_2^{\max (r_2,s_2)}...p_n^{\max (r_n,s_n)} \]

1.拓展歐幾里得演算法（exgcd）

洛谷P5656

用於求解形如\(ax+by=gcd(a,b)\)的不定方程特解。

當\(b=0\)時，可以看出\(gcd(a,b)=a\)，而

此時

\[\begin{cases} x=1 \\ y=0 \end{cases} \]

（實際上此時y大小不影響程式碼實現）

當\(b≠0\)時，遞迴求解\(exgcd(b,a\ mod\ b,x,y)\)，設

\[\begin{cases} a'=b \\ b'=a\% b \end{cases} \]

可以求得\(a'x'+b'y'=gcd(a,b)\)的一組特解，即\(x'\),\(y'\)。

所以得到了遞迴結束後\(x\)，\(y\)的表示式

\[\begin{cases} x=y' \\ y=x'-(a/b)*y' \end{cases} \]

證明如下：

程式碼：

void exgcd(int a,int b,int &x,int &y)
{
	if(!b)
	{
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,x,y);
	int p;
	p=x;
	x=y;
	y=p-(a/b)*y;//x=y',y=x'-(a/b)*y' 
	return;
}

還有一種更短的

void exgcd(int a,int b,int &x,int &y)
{
	if(!b)
	{
		x=1;y=0;
		return;
	}
	exgcd(b,a%b,y,x)//x=y',y=x'
	y-=(a/b)*x;//y=x'-(a/b)*y'
	return;
}

根據遞迴，我們可以知道這個x，y特解滿足

|x|+|y|最小

但是我們不滿足求這一組解

如果設

\(d = gcd(a,b)\)

那麼：

所以x，y的解可以寫成（x0+kb',y0-ka')

在此時，我們可以求出x，y最小非負整數解

分別是

xin=(x%b+b)%b;//最小非負整數解 

yin=(y%a+a)%a;//最小非負整數解

xin=x>0&&x%b!=0?x%b:x%b+b;//最小正整數解

yin=y>0&&y%a!=0?y%a:y%a+a;//最小正整數解

//最大整數解可以通過代入求出

當然，我們看到上面的求證過程中一直沒有出現用到
“\(ax+by\)右面是什麼”

那麼我們可以推廣：

設a,b,c為任意整數，若方程\(ax+by=c\)其中一個解是\((x0,y0)\)

則它的任意整數解可以寫成 \((x0+kb',y0-ka')\)

由此我們知道了任意整數解的求法，那\(ax+by=c\)的特解怎麼求呢？

這裡給出了一般性的做法，但為了編寫程式碼方便

我們一般這麼做

設

\[\begin{cases} g=gcd(a,b) \\ a'=a/g \\ b'=b/g \\ c'=c/g \\ \end{cases} \]

則

\[ax+by=c⇔a'x+b'y=c' \]

此時\(gcd(a',b')=1\),可以利用exgcd求出\(a'x'+b'y'=1\)的一組特解，繼而得出

\[\begin{cases} x0=c'x' \\ y0=c'y' \end{cases} \]

我們便求得了\(ax+by=c\)的一組特解。

這裡給出p5656的程式碼

//exgcd
#include <bits/stdc++.h>
#define int long long
using namespace std;
int gcd(int a,int b)
{
	if(b==0) return a;
	else return gcd(b,a%b);
}
void exgcd(int a,int b,int &x,int &y)
{
	if(!b)
	{
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,x,y);
	int p;
	p=x;
	x=y;
	y=p-(a/b)*y;//x=y',y=x'-(a/b)*y' 
	return;
}

signed main()
{
	
	int t;
	scanf("%lld",&t);
	while(t--)
	{
		int a,b,c,x=0,y=0,xin,xax,yin,yax,npa=0,g;//分別是x，y最小，最大正整數解 ,和正整數解的數量 
		scanf("%lld%lld%lld",&a,&b,&c);
		g=gcd(a,b);
		if(c%g!=0) printf("-1\n");//裴蜀定理 
		else
		{
			a/=g;
			b/=g;
			c/=g;//eg:求6x+15y=15:a：6/3=2,b:15/3=5,c:15/3=5
			//求2x'+5y'=1的一組解，x'=-2,y'=1
			//則原解為x'*c,x=-10,y=5;
			exgcd(a,b,x,y);//a'x+b'y=1
			x*=c;
			y*=c;
			//xin=(x%b+b)%b;最小非負整數解 
			xin=x>0&&x%b!=0?x%b:x%b+b;
			yax=(c-a*xin)/b;
			//yin=(y%a+a)%a;最小非負整數解 
			yin=y>0&&y%a!=0?y%a:y%a+a;
			xax=(c-b*yin)/a;
			if(xax>0)//yax>0也行
			npa=(xax-xin)/b+1;//正整數解數量
			//npa=(yax-yin)/a+1; 
			if(npa==0)
			{
				printf("%lld %lld\n",xin,yin);
			}
			else printf("%lld %lld %lld %lld %lld\n",npa,xin,yin,xax,yax);
            
		}
	}
	
	return 0;
	
	
}

2.線性同餘方程

對於形如 \(ax≡c(mod\ b)\) 的線性同餘方程，
根據模運算的定義，在方程左側新增一個\(by\)不會對結果造成影響，其實質就等價於\(ax+by=c\)的不定方程，利用exgcd求解便可。

注意：\(a≡c(mod\ b)\)有解的充要條件是：\(a-c\)是\(b\)的整數倍

例題：

洛谷P1082

轉換成\(ax\ mod\ b=1\)

轉換成移項可得\(ax+by=1\)（保證y是負數）

之後用exgcd求解

程式碼：

//轉化為求解ax+by=1 
#include <bits/stdc++.h>
using namespace std;
void exgcd(int a,int b,int &x,int &y)
{
	if(!b)
	{
		x=1;y=0;
		return;
	}
	exgcd(b,a%b,x,y);
	int temp=x;
	x=y;
	y=temp-(a/b)*y;
}
int main()
{
	int a,b,x,y;
	scanf("%d%d",&a,&b);
	exgcd(a,b,x,y);
	x=x>0&&x%b!=0?x%b:x%b+b; 
	printf("%d",x); 
	return 0;
} 

拓展性質：

$ax \equiv 1（mod\ b $稱為同餘方程的“逆”：a與b互質，且有唯一解

（注意：線性方程的唯一解是一組解）

它也是求解逆元的方法。。

3.乘法逆元

乘法逆元，一般用於求 \(\frac{a}{b} \pmod p\) 的值（\(p\)通常為質數），是解決模意義下分數數值的必要手段。

(1) exgcd

模數可以 不為質數

這個方法十分容易理解，而且對於單個查詢效率似乎也還不錯(尤其對於$ \bmod {p} $比較大的時候)。

這個就是利用拓歐求解 線性同餘方程$ ax \equiv c \pmod {b}\(的\)c=1\(的情況。我們就可以轉化為解\) ax + b*y = 1$這個方程。

求解這個方程的解。

而且這個做法還有個好處在於，當$ a \bot p$（互質），但 \(p\) 不是質數的時候也可以使用。

程式碼比較簡單：

void Exgcd(ll a, ll b, ll &x, ll &y) {
    if (!b) x = 1, y = 0;
    else Exgcd(b, a % b, y, x), y -= a / b * x;
}
int main() {
    ll x, y;
    Exgcd (a, p, x, y);
    x = (x % p + p) % p;
    printf ("%d\n", x); //x是a在mod p下的逆元
}

(2)費馬小定理

只適用於模數為質數的情況

若\(p\)為素數，\(a\)為正整數，且\(a\)、\(p\)互質。 則有\(a^{p-1} \equiv 1 (\bmod {p})\)。

另一個形式：

對於任意整數 $a $ ,有\(a^p \equiv \ a (mod \ p)\)

觀察第一個公式：

這個我們就可以發現它這個式子右邊剛好為 1 。

所以我們就可以放入原式，就可以得到：

\(a*x\equiv 1 \pmod p\)

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

$x \equiv a^{p-2} \pmod p $

所以我們可以用快速冪來算出 \(a^{p-2} \pmod p\)的值，這個數就是它的逆元了

ll fpm(ll x, ll power, ll mod) {
    x %= mod;
    ll ans = 1;
    for (; power; power >>= 1, (x *= x) %= mod)
        if(power & 1) (ans *= x) %= mod;
    return ans;
}

(3) 線性演算法

只適用於模數為質數的情況

用於求一連串數字對於一個\(\bmod p\)的逆元。

洛谷P3811

只能用這種方法，別的演算法都比這些要求一串要慢。

首先我們有一個,\(1^{-1}\equiv 1 \pmod p\)

然後設$ p=k*i+r,(1<r<i<p)$也就是 $ k$ 是 $ p / i$ 的商， $r $ 是餘數 。

再將這個式子放到\(\pmod p\)意義下就會得到：

\(k*i+r \equiv 0 \pmod p*\)

然後乘上\(i^{-1},r^{-1}\)就可以得到:

\(k*r^{-1}+i^{-1}\equiv 0 \pmod p\)

\(i^{-1}\equiv -k*r^{-1} \pmod p\)

\(i^{-1}\equiv -\lfloor \frac{p}{i} \rfloor*(p \bmod i)^{-1} \pmod p\)

於是，我們就可以從前面推出當前的逆元了。

注意：$ i ^{-1} * i ^{1} \equiv 1 $

#include <bits/stdc++.h>
#define ll long long
#define N 3000010
using namespace std;
ll inv[N];
int main()
{
	int n,p;
	scanf("%d%d",&n,&p); 
	inv[1]=1;
	printf("1\n");
	for(int i=2;i<=n;i++)
	{
		inv[i]=(ll)(p-p/i)*inv[p%i]%p;
		printf("%lld\n",inv[i]);
	}
	return 0;
 } 

（4） 階乘逆元法

只適用於模數為質數的情況

設$ f(i)=inv(i!)$, $ g(i)=i!\ $

則：$ f(i-1) = f(i) \times i $

• 證明:

$f(i-1)=\frac{1}{\ (i-1)\ !}=\frac{1}{i\ !}\times i =f(i)\times i $

假設要求 \([1,n]\) 中所有數的逆元

先求得 \([1,n]\) 中所有數的階乘

再用費馬小定理 求得\(f(n)\)的值

之後遞推出 \(f(1 \sim n)\) 的值

但是 \(inv(1! \sim n! )\) 並不是我們想要的答案

需要繼續轉化。

可知 : $inv(i) = inv(i!) \times(i-1)\ ! $

• 證明 :

  \(inv(i)=\frac{1}{i}=\frac{1}{i\ !}\times (i-1)\ ! = inv(i!)\times (i-1)!\)

按照上述方法轉換,
可得:

$ inv(i)=f(i)\times (i-1)!$

即得答案 。

#include<cstdio>
#define ll long long
using namespace std;
ll mul(ll a,ll b,ll mod) //快速冪模板
{
  ll ans=1;
  while(b)
    {
    	if(b&1) ans=ans*a%mod;
    	a=(a*a)%mod;
    	b>>=1;
    }
  return ans%mod;
}
ll n,p;
ll c[5000010]={1};
ll f[5000010];
int main()
{
  scanf("%lld%lld",&n,&p);
  for(int i=1;i<=n;i++)
    c[i]=(c[i-1]*i)%p;
    
  f[n]=mul(c[n],p-2,p); //獲得inv(n!)
  
  for(int i=n-1;i>=1;i--) //遞推階乘的逆元
    f[i]=(f[i+1]*(i+1))%p;
    
  for(int j=1;j<=n;j++) //轉化並輸出
    printf("%lld\n",(f[j]*c[j-1])%p);
}

4.素數篩法

（1）Eratosthenes 篩法 （埃拉託斯特尼篩法）

時間複雜度是\(O(nloglogn)\) 。

int Eratosthenes(int n) 
{
  int p = 0;
  for (int i = 0; i <= n; ++i) is_prime[i] = 1;
  is_prime[0] = is_prime[1] = 0;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime[p++] = i;  // prime[p]是i,後置自增運算代表當前素數數量
      if ((long long)i * i <= n)
        for (int j = i * i; j <= n; j += i)
          // 因為從 2 到 i - 1 的倍數我們之前篩過了，這裡直接從 i
          // 的倍數開始，提高了執行速度
          is_prime[j] = 0;  // 是i的倍數的均不是素數
    }
  }
  return p;
}

(2) Euler 篩法 （尤拉篩法）

時間複雜度\(O(n)\)。

洛谷P3383 線性篩模板

程式碼中，外層列舉 \(i = 1 \to n\)。對於一個 \(i\) ，經過前面的腥風血雨，如果它還沒有被篩掉，就加到質數陣列 \(Prime[]\) 中。下一步，是用 \(i\) 來篩掉一波數。

內層從小到大列舉\(Prime[j]\)。 \(i×Prime[j]\) 是嘗試篩掉的某個合數，其中，我們期望 \(Prime[j]\) 是這個合數的最小質因數 (這是線性複雜度的條件，下面叫做“篩條件”)。它是怎麼得到保證的？

\(j\) 的迴圈中，有一句就做到了這一點：

if(i % Prime[j] == 0)
	break; 

• 下面用 \(s(smaller)\) 表示小於 \(j\)的數，\(L(larger)\) 表示大於 \(j\) 的數。

• ① \(i\) 的最小質因數肯定是 \(Prime[j]\)。

（如果 \(i\) 的最小質因數是 \(Prime[s]\) ，那麼 \(Prime[s]\) 更早被列舉到（因為我們從小到大列舉質數），當時就要break）

既然 \(i\) 的最小質因數是 \(Prime[j]\)，那麼 \(i×Prime[j]\) 的最小質因數也是 \(Prime[j]\)。所以，\(j\) 本身是符合“篩條件”的。

• ② \(i×Prime[s]\) 的最小質因數確實是 \(Prime[s]\)。

（如果是它的最小質因數是更小的質數 \(Prime[t]\)，那麼當然 \(Prime[t]\) 更早被列舉到，當時就要break）

這說明 \(j\) 之前（用 \(i×Prime[s]\) 的方式去篩合數，使用的是最小質因數）都符合“篩條件”。

• ③ \(i×Prime[L]\) 的最小質因數一定是 \(Prime[j]\)。

（因為 \(i\) 的最小質因數是 \(Prime[j]\)，所以 \(i×Prime[L]\) 也含有 \(Prime[j]\) 這個因數（這是 \(i\) 的功勞），所以其最小質因數也是 \(Prime[j]\)（新的質因數 \(Prime[L]\) 太大了））

這說明，如果 \(j\) 繼續遞增（將以 \(i×Prime[L]\) 的方式去篩合數，沒有使用最小質因數），是不符合“篩條件”的。

#include <bits/stdc++.h>
using namespace std;
bool isprime[100000010];
int prime[6000010];
int cnt = 0;
void getprime(int n)
{
	memset(isprime,1,sizeof(isprime));
	isprime[1] = 0;//1不是素數
	for(int i=2;i<=n;i++)
	{
		if(isprime[i]) prime[++cnt] = i;
		for(int j=1;j<=cnt;j++)
		{
			if(i*prime[j]>n) break;
			isprime[i*prime[j]] = 0;
			if(i % prime [j] == 0) break;
		}
	}
}
int main()
{
	int n,q;
	scanf("%d%d",&n,&q);
	getprime(n);
	while(q--)
	{
		int k;
		scanf("%d",&k);
		printf("%d\n",prime[k]);
	}
	return 0;
}
	

5.擴充套件中國剩餘定理（exCRT）

洛谷P4777

給定\(n\)組非負整數\(a_i,b_i\)，求解關於\(x\)的方程組的最小非負整數解。

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

讓我們來改變一下格式：

\[\begin{cases} x+y_1a_1=b_1(1)\\ x-y_2a_2=b_2(2)\\ x-y_3a_3=b_3(3)\\ ... \end{cases} \]

把（1）（2）相減得：

\[y_1a_1+y_2a_2=b_1-b_2 \]

用\(\operatorname{exgcd}\)求解，不能解就無解。
然後我們可以解出一個最小正整數解\(y_1\)，帶入（1）得到\(x\)其中一個解：

\[x_0=b_1-a_1*y_1 \]

由於我們知道，\(y_1\)的全解，

\[y_1 '=y_1 + k*\frac{a_2}{\operatorname{gcd}(a_1,a_2)} \]

那麼x的全解是

\[x=b_1-a_1*y_1' \]

\[x=b_1-a_1*(y_1 + k*\frac{a_2}{\operatorname{gcd}(a_1,a_2)}) \]

\[x=b_1-a_1*y_1 - k*\frac{a_1*a_2}{\operatorname{gcd}(a_1,a_2)} \]

\[x=x_0+k\operatorname{lcm}(a_1,a_2) \]

由\(y_1\)的全解可導

即：\(x\equiv x_0(\ mod\ \operatorname{lcm}(a_1,a_2))\)

則：\(x+y_3\operatorname{lcm}(a_1,a_2)=x_0(4)\)

把（3）（4）再聯立

即可求解

6.關於取模問題

仍要記得開\(\operatorname{long long}\)！！

取模銘記“能取就取”！

公式：

\[(a+b)\ mod\ m=(a\ mod\ m+b\ mod\ m)\ mod\ m \]

\[(a-b)\ mod\ m = (a\ mod\ m-b\ mod\ m+m)\ mod\ m \]

\[(ab)\ mod\ m=\left[ (a\ mod\ m)\times (b\ mod\ m)\right]\ mod\ m \]

一定要記住乘法仍有可能爆\(\operatorname{int}\)

7.尤拉函式與尤拉定理

\(\phi (n)\)指不超過n並且與n互素的正整數的個數。

定理1：對於\(n=p_1^{a_1}p_2^{a_2}...p_n^{a_n}\)，有\(\phi (n)=\phi (p_1^{a_1})\phi (p_2^{a_2})...\phi (p_n^{a_n})\)

定理2：\(p\)為素數，則\(\phi (p)=p-1\).該定理充要。

定理3：\(p\)為素數，\(a\)是正整數，則\(\phi (p^a)=p^a-p^{a-1}\)

定理4：\(m,n\)為互質，則\(\phi (mn)=\phi (m)\phi (n)\).

定理5：\(p\)為奇數，則\(\phi (2n)=\phi (n)\).

定理6：\(n\)為大於2正整數，則\(\phi (n)\)是偶數.

定理7：\(n\)為正整數，則\(\sum _{d\mid n} \phi(d)=n\).

尤拉定理

對於任何兩個互質的正整數($\ a \bot m\ \() ，且\)a,m(m\geq 2)$

有\(a^{\phi(m)}\equiv 1(\mod m)\)

所以 \(a^b\equiv a^{b\bmod \varphi(m)}\pmod m\)

費馬小定理

尤拉定理中\(m\) 為質數時，\(a^{m-1}\equiv 1(\mod m)\)【尤拉定理+定理2】

應用：利用尤拉函式求不超過n且與n互素的正整數的個數，其次可以利用尤拉定理與費馬小定理來求得一個逆元，尤拉定理中的m適用任何正整數，而費馬小定理只能要求m是質數。

拓展尤拉定理

擴充套件尤拉定理無需 \(a,m\) 互質。

當 \(a,m\in \mathbb{Z}\) 時有：

\[a^b\equiv\begin{cases} a^b&,b<\varphi(m)\\ a^{b\bmod\varphi(m)+\varphi(m)}&,b\ge\varphi(m) \end{cases} \]