1. 数论基础
2. 线性同余方程的求解
3. 高次同余方程的求解
4. 同余方程组的求解
5. 组合数取摸
6. 幂取摸

数论基础

整除及其性质

def

对于式子 b = aq+r(a, b, q, r \in \mathbb{Z}, 0 \leq r \leq |a|)

当 r = 0 时，称 a 整除 b，记作 a | b

con

a|c, b | c, (a, b) = 1 \Rightarrow ab | c \\ a|bc, (a, b) = 1 \Rightarrow a | c \\ p | ab \Rightarrow p | a \ 或者 \ p | b

同余

def

对于式子 b = aq+r(a, b, q, r \in \mathbb{Z}, 0 \leq r \leq |a|)

则有

a \equiv r\mod b

con

a \equiv r \mod b \Rightarrow a | b - r \\ a \equiv b \mod n, a \equiv b \mod m \Rightarrow a \equiv b \mod [n, m] \\ d = (k, m), ka \equiv ka' \mod m \Rightarrow a \equiv a' \mod \frac{m}{d}

线性同余方程的求解

gcd 算法

算法作用

求 (a, b)

算法原理

欧几里得定理

定理内容

(a, b) = (b, a \mod b) (a \geq b)

定理证明

先证 (a, b) = (b, a - b)

设 d = (a, b)，dk_1 = a, dk_2 = b

则 a - b = (k_1 - k_2)d

\therefore d | (a-b)

\therefore d = (b, a - b)

可以发现欧几里得定理与上述定理本质相同

算法实现

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

当 b = 0 时，明显 a 就是 (a, b)

每次取摸，时间复杂度 O(\log a)

exgcd 算法

算法作用

判断与求不定方程

ax+by= (a, b)

算法引理

裴蜀定理

定理内容

ax+by=c

有整数解，当且仅当 (a, b) | c

定理证明

设 d = (a, b), dk_1 = a, dk_2 = b

则显然有

ax+by = dk_{1}x+dk_{2}y = d(k_1x + k_2y)

则

d | ax+by

则若原式有解，则必须满足 b|c

定理推论

(a, b) = 1 \Leftrightarrow ax+by=1 有整数解

算法原理

我们要求 ax + by = (a, b)

我们猜想，既然求解 (a, b) 是一个递归的过程，那么求解这个方程也可能会和递归有关，而且和 (a, b) 的递归基一样，当 b = 0 是，显然有 x = 1, y = 0

那么我们设

ax' + (a \text{ mod } b)y' = (b, a \text{ mod } b)

由欧几里得定理可以知道

(b, a \text{ mod } b) = (a, b)

所以

bx' + (a \text{ mod } b)y' = ax + by

注意到 a \mod b = a - \Big\lfloor\frac{a}{b}\Big\rfloor b

所以上式化为

bx' + (a - \Big \lfloor \frac{a}{b}\Big \rfloor b)y' = ax + by

整理得

b(x' - \Big\lfloor\frac{a}{b}\Big\rfloor y') + ay'

那么可以发现一组特解

x = y',\ y = x' - \Big\lfloor\frac{a}{b}\Big\rfloor y'

所以我们可以在递归右边的同时递归左边，解出这个方程

算法实现

pii exgcd(int a, int b) {
    if (b == 0) return mp(1, 0);
    pii tmp = exgcd(b, a % b);
    pii res = mp(tmp.second, tmp.first - a / b * tmp.second);
    return res;
}

时间复杂度 O(\log a)

高次同余方程的求解

BSGS

算法作用

求解高次方程

a^x \equiv b \mod m

其中

(a, m) = 1

算法原理

令 d = \sqrt m

设 i, j，使得 x = i * d - j\ (0 \leq i \leq m, 0 < j < m)

方程转化为

a^{i * d - j} \equiv b \mod m

打开得

a^{i*d} \equiv b*a^j \mod m

我们枚举 j = 0 \rightarrow d - 1，将 b\ast a^{j} 放入 \text{hash} 表中（键为 b\ast a^{j}，值为 i）

再枚举 i = 0 \rightarrow d，检查 \text{hash} 表中是否存在 a^{ i*d}，若有则 i, j 即为构成 x 的数，即可求解

感性的来说，就是分块思想在指数上的应用

算法实现

ll qpower(ll a, ll x, ll m) {
    ll res = 1;
    if (x == 0) return 1;
    while (x) {
        if (x & 1) res = res * a % m;
        a = a * a % m, x >>= 1;
    }
    return res% m;
}

ll bsgs(ll a, ll b, ll m) {
    unordered_map<ll, int> lst; lst.clear();
    ll d = ceil(sqrt(m));
    for (int j = 0; j <= d - 1; ++j) {
        lst[b * qpower(a, j, m) % m] = j;
    }
    ll tmp = qpower(a, d, m);
    if (tmp == 0) return b == 0 ? 1 : -1;
    for (int i = 0; i <= d; ++i) {
        ll val = qpower(tmp, i, m);
        if (lst.find(val) != lst.end() && i * d - lst[val] >= 0) return i * d - lst[val];
    }
    return -1;
}

时间复杂度 O(\sqrt m\log m)

同余方程组的求解

CRT

算法作用

求解同余方程组

\begin{cases} x \equiv a_1 \mod m_1 \\ x \equiv a_2 \mod m_2 \\ \vdots \\ x \equiv a_n \mod m_n \end{cases}

其中 \forall i, j(i \not = j) 都有 (m_i, m_j) = 1

算法原理

令

M = \prod_{i = 1}^n m_i, \ M_i = \frac{M}{m_i}, t_i = M_i^{-1} \mod m_i

则解为

x = \sum_{i = 1}^n a_iM_it_i

本质上是对于互不影响的 n 个方程进行依次回带的过程

算法实现

pii exgcd(int a, int b) {
    if (b == 0) return mp(1, 0);
    pii tmp = exgcd(b, a % b);
    pii res = mp(tmp.second, tmp.first - a / b * tmp.second);
    return res;
}

void CRT() {
    for (int i = 1; i <= n; ++i) {
        int inv = exgcd(prod / m[i], m[i]), inv = (inv % m[i] + m[i]) % m[i];
        ans += (((a[i] * prod) % m[i] + m[i]) % m[i] / m[i] * inv % m[i] + m[i]) % m[i];
    }
}

时间复杂度 O(n\log m)

exCRT

算法作用

求解同余方程组

\begin{cases} x \equiv a_1 \mod m_1\\ x \equiv a_2 \mod m_2\\ \vdots \\ x \equiv a_n \mod m_n \end{cases}

其中 m 没有任何特殊限制

算法原理

与 CRT 类似，exCRT 本质上仍然是一个回带的过程

对于两个方程

假设我们已经解出了 k 个形如

x \equiv a_i \mod m_i

的方程

且令

m = \text{lcm}_{i}^k m_i

我们有理由相信解出的 x 可以由

x_0 + t * m

的形式表达出

那么在考虑第 k + 1 个式子是，发现

x_0 + tm \equiv a_{k + 1} \mod m_{k + 1}

整理得

tm \equiv a_{k + 1} - x_0 \mod m_{k + 1}

即可由 exgcd 求解 t 值，然后更新 x_0 = x_0 + tm, 更新 m = m * m_{k + 1} / (m, m_{k + 1})

算法实现

#include <bits/stdc++.h>
#define pbb pair<__int128, __int128>
#define mp make_pair
using namespace std;
typedef __int128 bint;
const bint maxn = 2e5 + 10;
bint n, a[maxn], b[maxn], ans, prod;

bint read() {
    char ch = getchar(); bint x = 0, f = 1;
    while (ch < '0' || '9' < ch) { if (ch == '-') f = -f; ch = getchar(); }
    while ('0' <= ch && ch <= '9') { x = x * 10 + ch - '0', ch = getchar(); }
    return x * f;
}

void out(bint x) {
    if (x < 0) { putchar('-'); x = -x; }
    if (x > 9) out(x / 10);
    putchar(x % 10 + '0');
}

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

pbb exgcd(bint a, bint b) {
    if (b == 0) return mp(1, 0);
    pbb tmp = exgcd(b, a % b);
    pbb ans = mp(tmp.second, tmp.first - a / b * tmp.second);
    return ans;
}

bint mul(bint a, bint b, bint mod) {
    bint res = 0;
    while(b) {
        if(b & 1) res = (res + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return res;
}

void exCRT(bint k) {
    bint B = ((a[k] - ans) % b[k] + b[k]) % b[k];
    bint GCD = gcd(prod, b[k]);
    pbb sol = exgcd(prod, b[k]);
    sol.first = mul(sol.first, B / GCD, b[k]);
    ans += prod * sol.first, prod *= b[k] / GCD;
    ans = (ans + prod) % prod;
}

int main() {
    n = read();
    for (int i = 1; i <= n; ++i) b[i] = read(), a[i] = read();
    ans = a[1], prod = b[1];
    for (int i = 2; i <= n; ++i) {
        exCRT(i);
    }
    out(ans);
    return 0;
}

时间复杂度 O(n \log m)

组合数取模

卢卡斯定理

算法作用

求

\binom{n}{m}\ \text{mod}\ p

算法原理

经过一系列推到可知

\binom{n}{m} \equiv \binom{\big\lfloor\frac{n}{p}\big\rfloor}{\big\lfloor\frac{m}{p}\big\rfloor} * \binom{n\ \text{mod}\ p}{m\ \text{mod}\ p}\mod \ p

算法实现

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int maxn = 1e5 + 10;
ll n, m, p, prod[maxn], inv[maxn];

void clear() { memset(prod, 0, sizeof prod); memset(inv, 0, sizeof inv); }

ll lucas(ll n, ll m) {
    if (m < n) return 0;
    if (m < p) return prod[m] * inv[n] % p * inv[m - n] % p;
    return lucas(n % p, m % p) * lucas(n / p, m / p) % p;
}

void solve() {
    scanf("%lld %lld %lld", &n, &m, &p); clear();
    prod[0] = inv[0] = prod[1] = inv[1] =  1;
    for (int i = 2; i <= min(m + n, p); ++i) prod[i] = prod[i - 1] * i % p;
    for (int i = 2; i <= min(m + n, p); ++i) inv[i] = (p - p / i) * inv[p % i] % p;
    for (int i = 2; i <= min(m + n, p); ++i) inv[i] = inv[i - 1] * inv[i] % p;
    printf("%lld\n", lucas(m, n + m));
}

int main() { ll t;
    scanf("%lld", &t);
    while (t --) solve();
    return 0;
}

时间复杂度 O(p \log_{p}n)，\Omega(p)

扩展卢卡斯定理

\color{red}\text{待补}

幂取摸

欧拉定理

定理内容

a^{\varphi(m)} \equiv 1 \mod m

定理证明

易知 \{a, a^2, a^3 \cdots a^k\}(a^k \equiv 1 \mod p) 为 m 的简化剩余系的子群

又由于 m 的简化剩余系的阶为 \varphi(m)，且根据拉格朗日定理可知，k | \varphi(m)，设 s = \frac{\varphi(m)}{k}

所以 a^{\varphi(m)} \equiv a^{ks} \equiv (a^k)^s \equiv 1 ^ s \equiv 1 \mod p

所以得证

扩展欧拉定理

定理内容

a^b \equiv a^{(b\ \text{mod}\ \varphi(m)) + \varphi(m)} \mod m