bool isPrime(int a)
{
if (a == 1) return false;
for (int i = 2; (long long)i*i <= a; i++) if (a % i == 0) return false;
return true;
}
随着 n 的增大，试根法的效率会迅速下降。可以根据费马小定理：对于质数 p 和整数 a，有 a^{p-1}\equiv 1(\bmod\ {p})。如果 n 为质数，随机选择 a，则有 a^{n-1} \bmod n=1。再根据二次探测定理：若 p 为奇素数，则 x^2\equiv 1(\bmod\ p) 的解只能为 x\equiv 1(\bmod\ p) 和 x\equiv -1(\bmod\ p)。进行以下操作：
将 n-1 除以 k 次 2，变为为奇数 d，便于使用二次探测定理。
生成一个在 [1,n-1] 中的随机整数 a，计算 x = a^d\bmod n。
如果 x=1 或 x=n-1（即 x = -1），则继续判断。否则进行 k 次操作。每次操作将 x 更改为 x^2\bmod n，如果此时 x=n-1，则 n 为质数。
但是这种算法是一种概率性算法，存在误判的可能性。可以增加判断次数以减少误判的可能。时间复杂度为 O(k\log^3n)，其中 k 为判断的次数。
//快速幂
long long qpow(long long a, long long b, long long p)
{
long long res = 1;
a %= p;
while (b > 0)
{
if (b & 1) res = (res*a) % p;
a = (a*a) % p;
b >>= 1;
}
return res;
}
bool check(long long d, long long n)
{
long long a = 2 + rand()%(n-4); //取 2~n-1 中的随机数
long long x = qpow(a, d, n);
if (x == 1 || x == n - 1) return true;
while (d != n - 1)
{
x = (x * x) % n;
d *= 2;
if (x == 1) return false;
if (x == n - 1) return true;
}
return false;
}
bool is_prime(long long n, int k)
{
if (n == 1) return 0;
long long d = n - 1;
while (d % 2 == 0) d /= 2;
while (k--) if (!check(d, n)) return false; //测试 k 次
return true;
}
bool isp[N+5];
vector<int> prs;
void Eul_initp()
{
for (int i = 2; i <= N; i++) isp[i] = 1;
for (int i = 2; i <= N; i++)
{
if (isp[i]) prs.push_back(i);
for (int prm : prs)
{
if (i*prm > N) break;
isp[i*prm] = false;
if (i % prm == 0) break; //此时 i 的倍数一定能整出 prm，已经被标记过一定不为质数，可直接跳出循环
}
}
}
分解质因数
整数唯一分解定理：任何一个大于1的整数n都可以分解成若干个素因数的连乘积。
将一个正整数 x，遍历 [2,\lfloor\sqrt{x}\rfloor]，每次就可以找到 x 的 2 个因数，进行除法后，可分解质因数。由于从小至大遍历，如果 x 能被 i 整除，那么它一定能被 i 的最小质因子整除，而 x 在遍历到最小质因子时已经将这个因数除去，所以答案中不会存在合数。时间复杂度为 O(\sqrt{N})。
vector<int> facr(int x)
{
vector<int> res;
for (int i = 2; i * i <= x; i++)
{
while (x % i == 0)
{
x /= i;
res.push_back(i);
}
}
if (x != 1) res.push_back(x);
return res;
}
约数和倍数
正整数 x 的约数是指能整除 x 的数，倍数是指能被 x 整除的数。
最大公约数
```cpp
long long gcd(long long a, long long b)
{
while (b)
{
long long tmp = a;
a = b;
b = tmp % b;
}
return a;
}
```
当然，也可以使用 C++ <numeric> 中的 `gcd(a,b)` 函数。
## 最小公倍数
$a$ 和 $b$ 的最小公倍数 $\operatorname{lcm}(a,b)$ 为 $a$ 和 $b$ 共有的倍数中最小的数。设两个数的最大公约数为 $x$，$a=k_1x$，$b=k_2x$，则：
$$
\operatorname{lcm}(a,b)=k_1k_2x=\frac{a}{x}\frac{b}{x}x=\frac{ab}{x}=\frac{ab}{\gcd(a,b)}
$$
只要求出两数的最大公约数，即可求出最小公倍数。
```cpp
long long lcm(long long a, long long b) {return a*b/gcd(a,b);}
```
同样地，也可以使用 `lcm(a,b)` 函数。
# 快速乘与快速幂
在运用接下来的数论定理时，计算大数的乘法取模和幂的取模时，直接乘或使用 `pow` 函数后再取模会超出范围溢出，爆精度，所以使用快速乘和快速幂可以在计算过程中进行取模，从而避免数值溢出。
## 快速乘
我们可以不断将 $(a\times b) \bmod p$ 转为以下形式：
$$
\begin{cases}
(a\times\lfloor\frac{b}{2}\rfloor+a\times\lfloor\frac{b}{2}\rfloor+a)\bmod p & b\equiv1(\bmod\ 2) \\
(a\times\lfloor\frac{b}{2}\rfloor+a\times\lfloor\frac{b}{2}\rfloor)\bmod p & b\equiv0(\bmod\ 2) \\
\end{cases}
$$
直到 $b=0$ 时，返回结果 $0$。
以下为递归形式，时间复杂度为 $O(\log_2n)$，空间复杂度为递归时栈的深度，即 $O(\log_2n)$。
```cpp
long long fastmul(long long a, long long b, long long p)
{
if (b == 0) return 0;
long long ans = fastmul(a,b/2,p);
ans = (ans+ans)%p;
if (b & 1) return (ans+a)%p;
return ans;
}
```
如果想要降低空间复杂度，可以使用迭代的方法，使空间复杂度降至 $O(1)$。
```cpp
long long fastmul(long long a, long long b, long long p)
{
long long ans = 0;
a %= p;
while (b)
{
if (b & 1) ans = (ans+a)%p;
a = (a+a)%p;
b >>= 1;
}
return ans;
}
```
在代码中，`b & 1` 等同于 `b % 2 == 1`，`b >>= 1` 等同于 `b /= 2`，以上语句通过位运算提高运算效率。
## 快速幂
对于 $a^b\bmod p$，可以将它转换为以下形式：
$$
\begin{cases}
(a^{\lfloor\frac{b}{2}\rfloor}\times a^{\lfloor\frac{b}{2}\rfloor}\times a)\bmod p & b\equiv1(\bmod\ 2) \\
(a^{\lfloor\frac{b}{2}\rfloor}\times a^{\lfloor\frac{b}{2}\rfloor})\bmod p & b\equiv0(\bmod\ 2) \\
\end{cases}
$$
直到 $b=0$ 时，返回结果为 $1$。
以下递归形式的时间、空间复杂度为 $O(\log_2n)$。
```cpp
long long fastpow(long long a, long long b, long long p)
{
if (b == 0) return 1;
long long ans = fastpow(a,b/2,p);
ans = (ans*ans)%p;
if (b & 1) return (ans*a)%p;
return ans;
}
```
要想降低空间复杂度，可转为迭代形式。
```cpp
long long fastpow(long long a, long long b, long long p)
{
long long ans = 1;
while (b)
{
if (b & 1) ans = (ans*a)%p;
a = (a*a)%p;
b >>= 1;
}
return ans;
}
```
# 初等数论定理
## 欧拉定理
定义 $\phi(n)$ 为小于等于 $n$ 的正整数中与 $n$ 互质的数的个数，这被称为欧拉函数。若 $\gcd(a,p)=1$，则 $a^{\phi(p)}\equiv1(\bmod\ p)$。
例题：
求 $2025^{2025} \bmod 43$ 的值。
解：
上面数值过大，无法手动计算。
因为余数具有可乘性，所以 $2025^{2025}\equiv(2025\bmod43)^{2025}\equiv4^{2025}(\bmod\ 43)$。
因为 $4$ 与 $43$ 互质，此时可以使用欧拉定理：先手动计算小于等于 $43$ 且与 $43$ 互质的数的个数，结果为 $42$。
$4^{2025}=4^{48\times42+9}=(4^{42})^{48}\times2^9\equiv4^9\equiv64^3\equiv21^3\equiv16(\bmod\ 43)$。
所以最终结果为 $16$。
## 扩展欧拉定理
根据以上的例题可以发现以下规律：
$$
a^b\equiv
\begin{cases}
a^{b\bmod\phi(m)} & \gcd(a,m) = 1 \\
a^b & b < \phi(m) \\
a^{(b\bmod\phi(m))+\phi(m)} & \gcd(a,m)\neq1,b\ge\phi(m)
\end{cases}
(\bmod\ m)
$$
这就是扩展欧拉定理。
## 费马小定理
若 $p$ 为素数，且 $\gcd(a,p)=1$，则 $a^{p-1}\equiv1(\bmod\ p)$；即对于任意整数 $a$，有 $a^p\equiv a(\bmod\ p)$。 观察发现，费马小定理就是欧拉定理中 $p$ 为质数时的特殊情况。
应用：
- 素性测试：前面提到过。
- 求解同余方程（幂取余）：化简方程 $ax\equiv b(\bmod\ p)$（例：P2613 【模版】有理数取余、NEUOJ1391 Big Big Power）。
- 乘法逆元：见后面。
例题：求 $2^{2^n} \bmod 10007$ 的值（$1\le n\le10^6$）
解法：
根据费马小定理可知 $2^{10006}\equiv 1(\bmod\ 10007)$。由于余数具有可乘性和可加性，所以 $2^{2^n}\equiv1\times2^{2^n-2^n\bmod 10006}+2^{2^n\bmod 10006}\equiv2^{2^n\bmod 10006}(\bmod\ 10007)$，写两个快速幂求出 $k = 2^n \bmod 10006$ 的值和 $2^k \bmod 10007$ 的值。
## 裴蜀定理（贝祖定理）
若 $a$ 和 $b$ 不全为 0，对于任意整数 $x,y$，满足 $\gcd(a,b)|ax+by$，且存在整数 $x,y$ 使 $ax+by=\gcd(a,b)$。
逆定理：若 $a,b$ 不全为 0，正整数 $d$ 使它们的公因数，且存在整数 $x,y$ 使 $ax+by=b$，则 $d = \gcd(a,b)$。
应用：
- 扩展欧几里得
- 化简方程（例：P4549 【模板】裴蜀定理）
例题：
给定正整数 $a,b,m$，判断 $ax+by=m$ 是否有整数解。
解：如果 $\gcd(a,b)|ax+by$，即 $\gcd(a,b)|m$，则方程有整数解，否则没有。
## 威尔逊定理
对于一个大于 $1$ 的自然数 $n$，当且仅当 $n$ 为素数时，$(n-1)!\equiv-1(\bmod\ n)$。
逆定理：如果一个整数 $n>1$ 满足 $(n-1)!\equiv-1(\bmod\ n)$，那么 $n$ 为质数。
应用：
- 阶乘取模
- 素性测试
HDU 2973：
给定正整数 $n$，求以下表达式的值。
$$
S_n = \sum_{k=1}^{n}[\frac{(3k+6)!+1}{3k+7}-[\frac{(3k+6)!}{3k+7}]]
$$
解法：
如果 $3k+7$ 为质数，
$$
\begin{align}
A_k
&= [\frac{(3k+6)!+1}{3k+7}-[\frac{(3k+6)!}{3k+7}]] \\
&= [\frac{(3k+6)!+1}{3k+7}-(\frac{(3k+6)!+1}{3k+7}-1)] \\
&= 1
\end{align}
$$
如果 $3k+7$ 为合数，根据逆定理，则 $\frac{(3k+6)!+1}{3k+7}$ 与 $\frac{(3k+6)!}{3k+7}$ 整数部分相同，所以：
$$
\begin{align}
A_k
&= [\frac{(3k+6)!+1}{3k+7}-[\frac{(3k+6)!}{3k+7}]] \\
&= 0
\end{align}
$$
只需要枚举每一个 $k$，判断 $3k+7$ 是否为质数，进行对应的加法即可。
## 中国剩余定理
$$
\begin{cases}
x \equiv a_1(\bmod\ m_1) \\
x \equiv a_2(\bmod\ m_2) \\
\vdots \\
x \equiv a_n(\bmod\ m_n) \\
\end{cases}
$$
其中 $m_1,m_2,\ldots,m_n$ 两两互质。在模 $M=m_1m_2\ldots m_n$ 有解。具体来说，令 $M_i=M/m_i$，则解为 $\sum_{i=1}^{n}a_iM_iM_i^{-1} \bmod M$，其中 $M_i^{-1}$ 表示 $M_i$ 的乘法逆元（见后面）。
应用：
- 求同余方程组的解（例：P1495）
## 卢卡斯定理（Lucas 定理）
素数 $p$ 满足 $C_m^n \bmod\ p\equiv C_{m\bmod\ p}^{n\bmod\ p}\times C_{\lfloor m/p\rfloor}^{\lfloor n/p\rfloor}(\bmod\ p)$，其中 $C_{m}^{n}$ 表示组合数。
换一种形式：
对于非负整数 $m$ 和 $n$，以及素数 $p$，设 $m$ 和 $n$ 的 $p$ 进制展开分别为：
$$
m = m_k p^k + m_{k-1} p^{k-1} + \dots + m_0
$$
$$
n = n_k p^k + n_{k-1} p^{k-1} + \dots + n_0
$$
那么：
$$
C_m^n \equiv \prod_{i=0}^k C_{m_i}^{n_i} \pmod{p}
$$
例题：（P3807）给定正整数 $n,m,p$，求 $C_{n+m}^{n}\bmod p$ 的值。
解：
根据 Lucas 定理，$C_{n+m}^{n}\equiv C_{(m+n)\bmod\ p}^{n\bmod\ p}\times C_{\lfloor (m+n)/p\rfloor}^{\lfloor n/p\rfloor}(\bmod\ p)$，可以使用递归函数。
```cpp
long long Lucas(long long n, long long m, long long p)
{
if (m == 0) return 1;
return (C(n%p,m%p,p) * Lucas(n/p,m/p,p)) % p;
}
```
# 乘法逆元与同余方程
## 乘法逆元
如果 $ax\equiv1(\bmod\ b)$，则 $x$ 为 $a\bmod\ b$ 的逆元，记为 $a^{-1}$。只有 $a$ 与 $b$ 互质时，才存在 $x$。
### 求解方法 1：费马小定理
使用费马小定理的条件是 $b$ 为质数。一般在求解少量乘法逆元的情况下使用。
$$
a^{b-1}\equiv 1(\bmod\ b) \\
ax\equiv a^{b-1}(\bmod\ b) \\
x\equiv a^{b-2}(\bmod\ b)
$$
所以只需通过快速幂解出 $x$ 的值即可。时间复杂度为 $O(\log b)$。
```cpp
long long qpow(long long a, long long b, long long p)
{
long long ans = 1;
while (b)
{
if (b & 1) ans = (ans*a) % p;
a = (a*a) % p;
b >>= 1;
}
return ans;
}
int main()
{
int a,p;
cin >> a >> p;
cout << qpow(a,p-2,p) << endl;
return 0;
}
```
### 求解方法 2：扩展欧几里得算法
欧几里得算法即辗转相除法，用于求解两个数的最大公约数，重要公式为：$\gcd(a,b)=\gcd(a\bmod b,b)$。
扩展欧几里得算法用于求解 $ax+by=\gcd(a,b)$ 的解，在欧几里得算法中维护 $x$ 和 $y$ 的值。
当 $\gcd(a,b)=1$ 时，
$$
ax+by=1 \\
ax\bmod b+by\bmod b=1\bmod b \\
ax\equiv 1(\bmod\ b)
$$
所以就可将 $ax\equiv 1(\bmod\ b)$ 转为 $ax+by=1$ 的形式，并且只有 $a$ 与 $b$ 互质时，才可以使用扩展欧几里得算法。
设 $bx_0+(a\bmod b)y_0=\gcd(b,a\bmod b)=\gcd(a,b)$，则：
$$
bx_0+(a-b\lfloor\frac{a}{b}\rfloor)y_0=\gcd(a,b) \\
bx_0+ay_0-by_0\lfloor\frac{a}{b}\rfloor =\gcd(a,b) \\
ay_0+b(x_0-y_0\lfloor\frac{a}{b}\rfloor)=\gcd(a,b)
$$
对比 $ax+by=\gcd(a,b)$，可得到
$$
\begin{cases}
x = y_0 \\
y= x_0-y_0\lfloor\frac{a}{b}\rfloor
\end{cases}
$$
则本题可以递归 $x$ 和 $y$，直到除数 $b$ 为 $0$，此时 $ax=\gcd(a,0)=a$，可得 $x=1$，$y$ 可以为任意值，为了简化运算，通常返回 $y=0$。
时间复杂度为 $O(\log b)$。
```cpp
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a; //返回 a 和 b 的最大公约数 a
}
int x0,y0;
int t = exgcd(b,a%b,x0,y0);
x = y0;
y = x0-a/b*y0;
return t;
}
int get_inv(int a, int b)
{
int x,y;
exgcd(a,b,x,y);
return (x%b+b)%b; //将 x 转为正数并对 b 取余
}
int main()
{
int a,b;
cin >> a >> b;
cout << get_inv(a,b);
return 0;
}
```
### 求解方法 3：线性算法
如果需要大规模地计算乘法逆元，可以使用线性递推算法（例：P3811）。
如果 $a=1$，则 $x\equiv1(\bmod\ b)$，$x=a^{-1}$。
如果 $a>1$，设 $b=ka+r$（$0<r<a$），则：
$$
b\equiv0(\bmod\ b) \\
ka+r\equiv0(\bmod\ b) \\
r\equiv-ka(\bmod\ b) \\
r\times xr^{-1}\equiv-ka\times xr^{-1}(\bmod\ b) \\
x\equiv-kr^{-1}(\bmod\ b) \\
x\equiv-\lfloor\frac{b}{a}\rfloor\times(b\bmod a)^{-1}(\bmod\ b) \\
$$
由于这样操作会出现负数，可以修改为：
$$
x\equiv(b-\lfloor\frac{b}{a}\rfloor)\times(b\bmod a)^{-1}(\bmod\ b)
$$
最终时间复杂度为 $O(n)$。
```cpp
#include <iostream>
using namespace std;
long long inv[30000001];
int main()
{
int n;
long long b;
cin >> n >> b;
inv[1] = 1;
for (int i = 2; i <= n; i++) inv[i] = ((b-b/i)*inv[b%i]) % b;
for (int i = 1; i <= n; i++) cout << inv[i] << endl;
return 0;
}
```
## 同余方程
如果将乘法逆元中的 $1$ 改为其他数字，变为 $ax\equiv b(\bmod\ m)$，这称为线性同余方程。
当 $a$ 和 $m$ 互质时，
$$
ax\equiv b(\bmod\ m) \\
ax\times a^{-1}\equiv b\times a^{-1}(\bmod\ m) \\
x\equiv a^{-1}b(\bmod\ m)
$$
当 $a$ 和 $m$ 不互质时，根据同余方程的定理，方程有解当且仅当 $\gcd(a,m)|b$。设 $g = \gcd(a,m)$，将同余方程两边同时除以 $g$，变为：
$$
\frac{a}{g}x\equiv\frac{b}{g}(\bmod\ \frac{m}{g})
$$
此时 $\frac{a}{g}$ 与 $\frac{m}{g}$ 互质，可以转为 $a$ 和 $m$ 互质的情况。特解 $x_0\equiv(\frac{a}{g})^{-1}\frac{b}{g}(\bmod\ \frac{m}{g})$。
$$
(x+k\frac{m}{g})\equiv (\frac{a}{g})^{-1}\frac{b}{g}(\bmod\ \frac{m}{g}) \\
x\equiv x_0,x_0+\frac{m}{g},x_0+2\frac{m}{g},\ldots,x_0+(g-1)\frac{m}{g}(\bmod\ m)
$$
其中的乘法逆元可以使用扩展欧几里得算法求解。
因为 $k\frac{m}{g}$ 在模 $m$ 下的周期为 $g$，所以 $x$ 如果有解，那么解有 $g$ 个。
```cpp
#include <iostream>
#include <vector>
using namespace std;
int exgcd(int a, int b, int &x, int &y) //扩展欧几里得求最大公约数、乘法逆元
{
if (b == 0)
{
x = 1;
y = 0;
return a; //返回 a 和 b 的最大公约数 a
}
int x0,y0;
int t = exgcd(b,a%b,x0,y0);
x = y0;
y = x0-a/b*y0;
return t;
}
vector<int> solve(int a, int b, int m)
{
vector<int> ans;
int x,y;
int g = exgcd(a,m,x,y); //求 a,m 的最大公约数
if (b % g != 0) return ans;
int ap = a/g, bp = b/g, mp = m/g;
exgcd(ap,mp,x,y); //求 ap 模 mp 的乘法逆元
x = (x%mp+mp)%mp; //确保 x 为正数
int x0 = (bp*x)%mp; //求出特解
for (int k = 0; k < g; k++) ans.push_back((x0+k*mp)%m); //求出通解
return ans;
}
int main()
{
int a, b, m;
cin >> a >> b >> m;
vector<int> ans = solve(a,b,m);
for (int v : ans) cout << v << endl;
return 0;
}
```