NOI2019 不到一周了，开始总复习 颓废 咯 ~
一些数论算法以前已经写过了，本文主要用于总结一些杂七杂八的小算法，查缺补漏。
传送门 （说到传送门，Steam 夏日大促 Portal 2 只要 3 CNY 就是一玩起来就 3D 眩晕没法学习） ：
-
【知识总结】Miller-Robin和Pollard-Rho（咕 ~ ）
本文将涉及以下内容：
-
更相减损术、欧几里得算法及扩展欧几里得算法
-
逆元
-
中国剩余定理及扩展中国剩余定理
-
Lucas 定理
本文中所有字母如无特殊说明均默认为正整数。
更相减损术、欧几里得算法及扩展欧几里得算法
更相减损术
首先一个很显然的结论：
（证明略，反正看起来很显然就对了 —— 虽然这个显然的结论经常想不到 QAQ ）
于是可以就照着这个一直递归（或者迭代）下去，边界是 (gcd(a,0)=a) 。
欧几里得算法
这就是大家熟悉的 (gcd(a,b)=gcd(a mod b, b)) ，可以从上面那个结论推过来。同样可以通过递归的方式计算：
int gcd(const int a, const int b)
{
return b ? gcd(b, a % b) : a;
}
听说有时候更相减损术因为不用取模反而比这个还快？但是为什么我从来没见人写过更相减损术？？
扩展欧几里得算法
这玩意可不仅仅是算个最大公因数那么简单，它是用来解方程的 ……
问题：对于如下的方程（未知数是 (x) 和 (y) ），求任意一组整数解：
令 (m=a-bcdotlfloorfrac{a}{b} floor) （即 (a) 模 (b) ），设 (x') 和 (y') 满足以下方程：
因为 (gcd(a,b)=gcd(m,b)) ，所以：
把 (m) 按照定义展开：
因此 (x=y') ， (y=x'-lfloorfrac{a}{b} floor 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}) 即可。
另外，线性同余方程 (axequiv bpmod 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=1pmod m) ，则 (a) 是 (b) 在模 (m) 意义下的逆元。这里介绍它的三种求法：
费马小定理
对于 质数 (p) 和任意正整数 (a) ，有 (a^{p-1}=1pmod 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);
}
扩展欧几里得
相当于同余方程 (axequiv 1pmod 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) 的逆元。
边界是 (1^{-1}=1) 。
证明：记 (t=lfloorfrac{m}{a} floor) ，(k=(mmod a)) ，则 (at+k=0pmod m) 。
两边同时乘上 (a^{-1}cdot k^{-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) 。
对于如下方程组：
其中 (m_1) 到 (m_n) 两两互质。
我们可以直接构造出一个解：
其中 (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} mod m_j=0(i eq 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) 不保证两两互质了。
考虑如果只有两个方程如何求解：
对于第一个方程而言，解的形式是 (x=a_1+km_1) 。带入第二个方程，得到：
这个方程中只有 (k) 是未知的。用扩展欧几里得解出来一个 (k_0) ，则任意 (km_1=k_0m_1+zm_2) 都是合法的（实在不知道用什么字母就用「嘴子」的第一个字母 (z) 好了），所以 (zm_2=ucdotmathrm{lcm}(m_1,m_2)) (实在不知道用什么字母就用「嘴子」的第二个字母 (u) 好了）。现在，(x) 的表现形式为 (a_1+k_0m_1+ucdotmathrm{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 以后再补。
其中 (p) 是质数。