【知识总结】数论全家桶

淩亂°似流年 2021-11-12 00:38 390阅读 0赞

NOI2019 不到一周了，开始总复习 颓废 咯 ~

一些数论算法以前已经写过了，本文主要用于总结一些杂七杂八的小算法，查缺补漏。

传送门 （说到传送门，Steam 夏日大促 Portal 2 只要 3 CNY 就是一玩起来就 3D 眩晕没法学习） ：

*  [【知识总结】线性筛\_杜教筛\_Min25筛][Min25]
 *  [【知识总结】扩展卢卡斯定理（exLucas）][exLucas]
 *  【知识总结】Miller-Robin和Pollard-Rho（咕 ~ ）

## 本文将涉及以下内容： ##

*  更相减损术、欧几里得算法及扩展欧几里得算法
 *  逆元
 *  中国剩余定理及扩展中国剩余定理
 *  Lucas 定理

本文中所有字母如无特殊说明均默认为正整数。

## 更相减损术、欧几里得算法及扩展欧几里得算法 ##

### 更相减损术 ###

首先一个很显然的结论：

\\\[ \\gcd(a,b)=\\gcd(a,a-b)(a\\geq b) \\\]

（证明略，反正看起来很显然就对了 —— 虽然这个显然的结论经常想不到 QAQ ）

于是可以就照着这个一直递归（或者迭代）下去，边界是 \\(\\gcd(a,0)=a\\) 。

### 欧几里得算法 ###

这就是大家熟悉的 \\(\\gcd(a,b)=\\gcd(a \\bmod b, b)\\) ，可以从上面那个结论推过来。同样可以通过递归的方式计算：

int gcd(const int a, const int b) { return b ? gcd(b, a % b) : a; }

听说有时候更相减损术因为不用取模反而比这个还快？但是为什么我从来没见人写过更相减损术？？

### 扩展欧几里得算法 ###

这玩意可不仅仅是算个最大公因数那么简单，它是用来解方程的 ……

问题：对于如下的方程（未知数是 \\(x\\) 和 \\(y\\) ），求任意一组整数解：

\\\[ax+by=\\gcd(a,b)\\\]

令 \\(m=a-b\\cdot\\lfloor\\frac\{a\}\{b\}\\rfloor\\) （即 \\(a\\) 模 \\(b\\) ），设 \\(x'\\) 和 \\(y'\\) 满足以下方程：

\\\[bx'+my'=\\gcd(m, b)\\\]

因为 \\(\\gcd(a,b)=\\gcd(m,b)\\) ，所以：

\\\[bx'+my'=ax+by\\\]

把 \\(m\\) 按照定义展开：

\\\[\\begin\{aligned\} ax+by&=bx'+(a-b\\cdot\\lfloor\\frac\{a\}\{b\}\\rfloor)y'\\\\ &=bx'+ay'-b\\cdot\\lfloor\\frac\{a\}\{b\}\\rfloor y'\\\\ &=ay'+b(x'-\\lfloor\\frac\{a\}\{b\}\\rfloor y') \\end\{aligned\}\\\]

因此 \\(x=y'\\) ， \\(y=x'-\\lfloor\\frac\{a\}\{b\}\\rfloor y'\\) 。于是我们得到了一组解。

然而，这是一个不定方程，解并不是唯一的。对于其中一组解 \\(\\\{x\_0,y\_0\\\}\\) ，有通解 \\(\\\{x\_0+kb,y\_0-ka\\\}\\) 。

形如 \\(ax+by=c\\) 都可以通过转化来用扩展欧几里得解决。设 \\(g=\\gcd(a,b)\\) ，如果 \\(g\\) 不能整除 \\(c\\) ，则原方程无整数解；否则，先解出 \\(ax+by=g\\) 的解，然后给解乘上 \\(\\frac\{c\}\{g\}\\) 即可。

另外，线性同余方程 \\(ax\\equiv b\\pmod p\\) 也可用扩展欧几里得解决，转化成 \\(ax+py=b\\) 即可。

int exgcd(const int a, const int b, int &x, int &y) { if (!b) { x = 1, y = 0; return a; } int xx, yy, tmp = exgcd(b, a % b, xx, yy); x = yy; y = xx - a / b * yy; return tmp; }

## 逆元 ##

定义：若 \\(ab=1\\pmod m\\) ，则 \\(a\\) 是 \\(b\\) 在模 \\(m\\) 意义下的逆元。这里介绍它的三种求法：

### 费马小定理 ###

对于 **质数** \\(p\\) 和任意正整数 \\(a\\) ，有 \\(a^\{p-1\}=1\\pmod p\\) 。因此在模 \\(p\\) 意义下 \\(a\\) 的逆元是 \\(a^\{p-2\}\\) ，直接快速幂计算即可。

int power(int a, int b) { int ans = 1; while (b) { if (b & 1) ans = (ll)ans * a % p; a = (ll)a * a % p; b >>= 1; } return ans; } int inv(const int a, const int p) { return power(a, p - 2); }

### 扩展欧几里得 ###

相当于同余方程 \\(ax\\equiv 1\\pmod m\\) ，直接用扩展欧几里得解方程即可。

exgcd 的实现见上文。

int inv(const int a, const int p) { int x, y, tmp = exgcd(a, p, x, y); if (tmp > 1) return -1; return (x % p + p) % p; }

### 线性递推逆元 ###

（这玩意我学了 N 遍甚至拿这个式子当电脑桌面一个月都没记住，很气。）

记 \\(a^\{-1\}\\) 表示 \\(a\\) 的逆元。

\\\[ a^\{-1\}= m-\\lfloor\\frac\{m\}\{a\}\\rfloor\\cdot (m\\bmod a)^\{-1\}\\pmod m\\\]

边界是 \\(1^\{-1\}=1\\) 。

证明：记 \\(t=\\lfloor\\frac\{m\}\{a\}\\rfloor\\) ，\\(k=(m\\bmod a)\\) ，则 \\(at+k=0\\pmod m\\) 。

两边同时乘上 \\(a^\{-1\}\\cdot k^\{-1\}\\) ，得到：

\\\[ tk^\{-1\}+a^\{-1\}=0 \\\]

即：

\\\[a^\{-1\}=-tk^\{-1\}\\\]

void init() { inv[1] = 1; for (int i = 2; i < p; i++) inv[i] = (p - (ll)p / i * inv[p % i] % p) % p; }

## 中国剩余定理及扩展中国剩余定理 ##

### 中国剩余定理 ###

大名鼎鼎的 CRT (Chongqing Rail Transit) (Chinese Remainder Theorem) 。

对于如下方程组：

\\\[\\begin\{cases\} x=a\_1\\pmod \{m\_1\}\\\\ x=a\_2\\pmod \{m\_2\}\\\\ \\cdots\\\\ x=a\_n\\pmod \{m\_n\}\\\\ \\end\{cases\}\\\]

其中 \\(m\_1\\) 到 \\(m\_n\\) 两两互质。

我们可以直接构造出一个解：

\\\[x=\\sum\_\{i=1\}^\{n\}a\_i\\frac\{M\}\{m\_i\}\\cdot \\mathrm\{inv\}(\\frac\{M\}\{m\_i\},m\_i)\\\]

其中 \\(M=\\prod\_\{i=1\}^n m\_i\\) ，\\(\\mathrm\{inv\}(a,b)\\) 表示 \\(a\\) 在模 \\(b\\) 意义下的逆元。由于模数两两互质，此时逆元一定存在。

对于任意一个 \\(i\\) ，当模 \\(m\_i\\) 时 \\(\\frac\{M\}\{m\_i\}\\cdot \\mathrm\{inv\}(\\frac\{M\}\{m\_i\},m\_i)\\) 是 \\(1\\) （逆元的定义），否则为 \\(0\\) （\\(\\frac\{M\}\{m\_i\} \\bmod m\_j=0(i\\neq j)\\)），因此这个 \\(x\\) 是原方程组的解。可以证明在模 \\(M\\) 的意义下有且只有这一个解。

int CRT(const int n) { int M = 1, ans = 0; for (int i = 0; i < n; i++) M *= p[i]; for (int i = 0; i < n; i++) ans = (ans + (ll)a[i] * inv(M / p[i], p[i]) % M * (M / p[i]) % M) % M; return ans; }

### 扩展中国剩余定理 ###

仍然是上面的问题，只是 \\(m\_1\\) 到 \\(m\_n\\) 不保证两两互质了。

考虑如果只有两个方程如何求解：

\\\[\\begin\{cases\} x=a\_1\\pmod \{m\_1\}\\\\ x=a\_2\\pmod \{m\_2\} \\end\{cases\}\\\]

对于第一个方程而言，解的形式是 \\(x=a\_1+km\_1\\) 。带入第二个方程，得到：

\\\[a\_1+km\_1=a\_2\\pmod \{m\_2\}\\\]

这个方程中只有 \\(k\\) 是未知的。用扩展欧几里得解出来一个 \\(k\_0\\) ，则任意 \\(km\_1=k\_0m\_1+zm\_2\\) 都是合法的（实在不知道用什么字母就用「嘴子」的第一个字母 \\(z\\) 好了），所以 \\(zm\_2=u\\cdot\\mathrm\{lcm\}(m\_1,m\_2)\\) (实在不知道用什么字母就用「嘴子」的第二个字母 \\(u\\) 好了）。现在，\\(x\\) 的表现形式为 \\(a\_1+k\_0m\_1+u\\cdot\\mathrm\{lcm\}(m\_1,m\_2)\\) ，我们合并出了一个新的方程：

\\\[x=a\_1+k\_0m\_1\\pmod\{\\mathrm\{lcm\}(m\_1,m\_2)\}\\\]

你看着像不像风靡机房的小游戏 全民合车 全民漂移（划掉）。

于是这样合并 \\(n\\) 次就完了。

顺便这里有一个叫「龟速乘」的小 trick ，利用整数溢出是循环溢出的特性制成 ……

ll mul(ll a, ll b, const ll p) { a %= p, b %= p; return ((a * b - ll((long double)a * b / p) * p) % p + p) % p; } ll excrt() { ll M = m[1], x = a[1]; for (int i = 2; i <= n; i++) { ll tmp = lcm(M, m[i]), res = solve(M, a[i] - x, m[i]); if (res == -1) return -1; x = (x + mul(res, M, tmp)) % tmp; M = tmp; } return x; }

## Lucas 定理 ##

就是一个公式，证明什么的 NOI 以后再补。

\\\[C\_n^m=C\_\{\\lfloor\\frac\{n\}\{p\}\\rfloor\}^\{\\lfloor\\frac\{m\}\{p\}\\rfloor\}\\cdot C\_\{n\\bmod p\}^\{m\\bmod p\} \\pmod p\\\]

其中 \\(p\\) 是质数。

转载于:https://www.cnblogs.com/zyt1253679098/p/11179050.html

[Min25]: https://www.cnblogs.com/zyt1253679098/p/10428767.html
[exLucas]: https://www.cnblogs.com/zyt1253679098/p/9736741.html