最大公约数、扩展欧几里得算法、逆元
数论基础
数论是纯数学的一个研究分支，主要研究整数的性质。初等数论包括整除理论、同余理论、连分数理论。这一篇主要记录的是同余相关的基础知识。
取模
取模是一种运算，本质就是带余除法，运算结果就是余数。取模运算结果的符号由被模数（被除数）决定。
取模运算的性质
例题1
找规律，对于mod n来说，最多n²项就会出现重复，找出重复项即可得到循环周期，然后找对应的项即可。
例题2：hdu 1212
题意是给一个位数不超过1,000的正整数A和一个大小不超过100,000的正整数B，求A%B。
这里只需要用到上面的两条性质，把A截断，从高到低模拟做除法即可。
#include <iostream>
string s; int mod, ans;
int main() {
while(std::cin >> str >> mod) {
ans = 0;
for(int i = 0; s[i]; i++) ans = (ans * 10 + (s[i] - '0')) % mod;
std::cout << ans << std::endl;
}
return 0;
}
GCD
GCD即Greatest Common Divisor，最大公约数。最大公约数即能同时整除给定两个整数的最大正整数。特殊地：gcd（a，0）=a。求解最大公约数通常使用辗转相除法。下面是辗转相除法的代码实现：
typedef long long LL;
LL gcd(LL a, LL b) {
return b ? gcd(b, a % b) : a;
}
辗转相除法得到的余数序列增长速度比斐波那契数列更快，已知斐波那契的增长是指数级别，则辗转相除法的复杂度是对数级别。
辗转相除法的证明：
1.设两个数a、b（a＞b），他们的最大公约数为gcd（a，b），r = a % b，k = （a - a % b）/ b，那么证明辗转相除法，即是要证明：gcd（a，b）=gcd（b，r）。首先我们令c = gcd（a，b），那么肯定存在互质的整数m，n，使得a = mc，b = nc。那么有r = a - kb = mc - knc = (m - kn)c，根据这条式子，c也是r的因数。回过头看，如果能证明m-kn和n互质，那么就可以说明c=gcd（（m-kn）c，nc）=gcd（b，r），所以把问题再次转化为：求证m-kn和n互质。
2.反证法证明m-kn和n互质：假设m-kn和n不互质，用数学语言描述为：假设存在整数x，y，d，其中d＞1，使得m-kn=xd，n=yd。那么有m=kn+xd=kyd+xd=(ky+x)d，从而a=mc=(ky+x)cd，b=nc=ycd，则gcd（a，b）=cd≠c，与前设矛盾。故m-kn和n互质得证，也就证明了gcd（a，b）=gcd（b，r）。
唯一分解定理
对于任意大于2的正整数X，它总能写成如下形式：
即质因数分解。且这个分解唯一。
LCM
LCM即Least Common Multiple，最小公倍数。
LL lcm(LL a, LL b) {
return a / gcd(a, b) * b;
}
例题：hdu1788
题意是求最小的正整数N，满足除以每一个给定的数Mi，其结果均为Mi-a。作如下转化：
即N是Mi的倍数，那么题意也就是求所有M的最小公倍数，注意数字范围即可。
同余
对于两个不同的数a，b，如果有a % p = b % p（p＞1），那我们就称a和b模p同余，记作：
同余的性质
EXGCD
贝祖定理（Bezouts Identity）：若设a,b是不全为0的整数，则存在整数x,y，使得ax+by=gcd（a,b），(a,b)代表最大公因数，则设a,b是不全为零的整数，则肯定存在整数x,y，使得ax+by=(a，b)。这个式子称为贝祖等式。
EXtend GCD 即扩展欧几里得算法，它可以用于解关于x，y的贝祖等式。
EXGCD的可行性
当a=0时，ax+by=gcd（a，b）=gcd（0，b）=b，这时可以解出y=1而x为任意。同理可得b=0时，解得x=1而y为任意。当a，b都不为0时：
a或b为0即迭代求解的出口，迭代过程用代码表示如下：
typedef long long LL;
// ver 1
// 返回值为gcd（a，b），x和y即求出的特解
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(b, a % b, x, y);
//这里通过迭代求出了一组解x，y，这组解对应上面推导过程中的x2和y2。
LL temp = x;
x = y;
y = temp - a / b * y;
return ans;
}
// ver 2
// 在熟悉了推导过程之后，我们可以利用引用的性质得到更为简洁的ver2写法
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(b, a % b, y, x);
y -= a / b * x;
return ans;
}
EXGCD的应用
求关于x，y的方程ax+by=c的一组解
这里只需要对贝祖等式稍作变换即可。假设EXGCD求出的方程ax+by=gcd（a，b）的一组解为x1和y1。在方程两边同时乘上一个数m，使得c=m*gcd（a，b），这里也就要求c是gcd（a，b）的倍数，即c%gcd（a，b）=0。这也是方程有解的条件。此时方程为：
对于ax+by=gcd（a，b），其参数解为：
那么对比方程ax+by=c，可以得到其特解为：
那么其通解为：
对于任意一个确定的t，都有一组确定的解与之对应，只需要根据需要找出对应的解即可。
例如求满足ax+by=c的X的最小非负整数解，那么在得到X的通解之后只需调整t，调整过程用伪代码表示为：
S = b / gcd(a, b);
X = (x0 % S + S) % S;
// 由于x0可能是负数，故在模之后还要加上S在取一次模。
// 同理可得Y的最小非负整数解
T = a / gcd(a, b);
Y = (y0 % T + T) % T;
完整的exgcd求解方程ax+by=c的X和Y的最小非负整数解得代码如下：
LL exgcd(LL a, LL b, LL c, LL& x, LL& y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
LL ans = exgcd(a, b, c, y, x);
y -= a / b * x;
LL S = b / ans, T = a / ans;
x = (x % S + S) % S;
y = (y % T + T) % T;
return ans;
}
求关于x的方程ax≡b（mod m）的最小非负整数解
ax≡b（mod m）即ax%m=b%m，即求ax+my=b%m（取模其实就相当于减去（加上）了若干个模数）。但只有a和m互质时有唯一解，否则无解。转化为ax+my=b（假设b<m）后，套用exgcd即可。
求a关于模数p的乘法逆元
设x是a关于p的乘法逆元，那么有ax≡1（mod p），这是第二个应用的特殊情况，显然也可以通过类似方法求得。
逆元
Inverse Element，逆元，推广了加法中的加法逆元和乘法中的倒数。直观地说，它是一个可以取消另一给定元素运算的元素。a关于模p的逆元存在的条件是gcd（a，p）=1。
为什么需要逆元呢？
这里，我们把除法转化为乘法，就可以运用取模运算的性质：(a * b) % c = (a % c * b % c) % c，优化算法。
逆元的求法
EXGCD求法
给定模数m，求a的逆元相当于求解关于x的方程ax≡1(mod m)。这个方程可以转化为ax-my=1 。用EXGCD可以求得一组x，y和gcd。检查gcd是否为1（因为EXGCD是把问题转化为求解方程ax-my=gcd（a，m），显然只有gcd（a，m）=1时才能转化），gcd不为1则说明逆元不存在。在能够成功求解的情况下，把解调整x到0~m-1的范围中即可。
费马小定理
在模p为质数的情况下，有费马小定理：
那么:
然后用快速幂求出逆元即可。
拓展：快速幂
// a是底数，b是指数
LL pow(LL a, LL b, LL mod) {
LL ans = 1;
while(b) {
if(b&1) ans = ans * a % mod;
b >>= 1, a = a * a % mod;
}
return ans;
}
拓展：慢速乘，防止快速幂过程中乘法运算爆long long，只需把快速幂中的乘法替换为mul（）即可，一般情况用不上。
// a为因数1，b为因数2
LL mul(LL a, LL b, LL mod) {
LL ans = 0;
while(b) {
if(b&1) ans = (ans + a) % mod;
b >>= 1, a = (a + a) % mod;
}
return ans;
}
欧拉定理
费马小定理其实是欧拉定理的特殊情况。在模数p不为质数时，有：
同理可以得到：
线性逆元表
有时我们需要快速得出1~p-1的所有逆元，这个时候我们需要一种O（n）的方法，这就是线性逆元表。
首先，1关于任意模p的逆元的逆元都是1。设p = k * i + r（r < i ，1 < i < p），那么有：
利用逆元性质，两边同时乘上i的逆元以及r的逆元，得到：
移项得：
显然p%i是小于i的，那么我们就可以利用之前的结果在O（1）时间算出单独的逆元。
用代码表示即：
// inv[i]对应i的逆元
LL inv[maxn];
void CalInv() {
// 0没有逆元，故不初始化
inv[1] = 1;
for (int i = 2; i < maxn; i++)
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
利用
我们还可以在线性时间求出1~min(n，p)的阶乘的逆元，代码如下：
LL fac[maxn], inv[maxn];
void CalFacInv() {
// 0的阶乘是1
fac[0] = fac[1] = 1;
for (int i = 2; i < maxn; i++)
fac[i] = fac[i - 1] * i % mod;
inv[maxn - 1] = QPow(fac[maxn - 1], mod - 2);
// 注意边界是≥
for (int i = maxn - 2; i >= 0; i--)
inv[i] = inv[i + 1] * (i + 1) % mod;
}
例题：hdu1576
题意即把除法转化为逆元乘法。3种方法的代码：
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
LL n, B;
const LL mod = 9973;
LL exgcd(LL a, LL b, LL &x, LL &y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
LL ans = exgcd(b, a % b, y, x);
y -= a/b*x;
return ans;
}
LL QPow(LL x, LL n)
{
LL res = 1;
while(n)
{
if(n & 1) res = res * x % mod;
x = x * x % mod;
n >>= 1;
}
return res;
}
LL inv[10000];
void CalInv() {
inv[1] = 1;
for (LL i = 2; i < 10000; i++) {
inv[i] = (mod * mod - mod / i * inv[mod % i]) % mod;
}
}
int main() {
CalInv();
int t;
cin >> t;
while (t--) {
cin >> n >> B;
// exgcd
// LL x, y;
// exgcd(B, mod, x, y);
// cout << (x + mod) * n % mod << endl;
// QPow
// cout << QPow(B, mod - 2) * n % mod << endl;
// 逆元打表
cout << inv[B % mod] * n % mod << endl;
}
return 0;
}
