2024年12月一些算法比赛（二）

[[ slider_text ]]

[[ item.c ]]

2024年12月一些算法比赛（二）

发布时间 2024-12-21 作者：zhangminchen 来源：算法小站

并查集 组合数学 凸包 dp优化 数论逆元

## AtCoder Beginner Contest 383
### E
**[E - Sum of Max Matching](https://atcoder.jp/contests/abc383/tasks/abc383_e)**
给定一个$$N, M$$的无向图，边带权，定义一条路径的权值，是路径上 weight 最大的边的权重
定义$$f(x, y)$$为顶点$$x \to y$$路径上，最小边的权重
现在给你两个序列，$$A, B$$，长度都是$$K$$
其中$$A$$顺序不能变，但$$B$$顺序可以变，要求$$\displaystyle\sum\_{i=1}^{K}f(A\_i, B\_i)$$最小化

#### 算法分析
> 类似 kruskal 思想

想法并不难，对边从小到大排序，然后考虑每条边，$$(u, v)$$，如果他们不在一个集合中，并查集合并
合并的时候，考虑`findA, findB`，把`findA, findB`中需要的点（需要的点就是，在序列$$A, B$$中出现的那些点）都匹配了即可

> 难点是怎么实现？

因为是匹配，考虑$$(u, v)$$合并，`pA = find(u),  pB = find(v)`
一般的匹配问题，需要维护**两类 item 分别有多少个？**，这里也同样
每个集合维护**它能提供多少个 A 类数和 B 类数？**（集合`A`中出现的数称为`A 类数`）

因为题目中保证了$$A_i \neq B_j$$，一开始的时候，某个节点`u`，要么都不在`A, B`集合中
要么它只能是`A类数`，或者是`B类数`，形如`needA(u) += 1, needB(u) = 0`
这样对于每个点`u`，`needA(u) = ?`，表示它要作为`A`类数被匹配`?`次
而当前`needB(u) = 0`，表示当前集合只能完成`0`次匹配，这样需要和其他集合合并

合并的时候，`fa( pB ) = pA`，同时要令`now = pA`
`needA( now ) = needA( pA ) + needA( pB ),  needB( now ) = needB( pA ) + needB( pB )`
`needA, needB`也要同步合并

> 从匹配 + 点的角度入手，我们先预处理序列，把需要用到这个点的都算上，这个点要被用多少次？
先一次性全部加上，这样$$(u, v)$$并查集合并的时候，就把`A, B`序列中，所有$$u, v$$出现的次数都算完了
之后如果还未匹配的，接着做并查集合并，如果`find(u) != find(v)`，继续匹配

本题的算法是正确的，因为输入保证了$$A_i \neq B_j$$，并查集合并之后，要么`needA = 0`，要么`needB = 0`

### F
**[F - Diversity](https://atcoder.jp/contests/abc383/tasks/abc383_f)**
有$$n$$种商品，每一种$$p_i$$元，同时有价值$$u_i$$，颜色$$c_i$$，现在要从$$n$$种商品选一个子集
可能一个都不选，然后你要最大化$$S + T\cdot K$$，其中$$K$$给定，$$S$$是所选的商品的价值和，$$T$$是有多少个不同的颜色
当然，商品的价格总和不能超过$$X$$

#### 算法分析

**记念一下，这个问题自己独立想出了一个新做法**

> 不难发现，如果没有$$T$$的限制，问题就是一个`01背包`

朴素的 dp 可以写出来
```bash
for i = 1 to N
	for x = 0 to X, 剩余的钱数
		dp(i, x) <-- dp(i-1, x)
		dp(i, x) <-- dp(i-1, x-p(i)) + u(i) + ( 若 c(i) 出现过 ？0 : k )
```

> 注意到，这个问题中，每一个颜色，只会在它第一次被选的时候，影响 dp 值，加上$$K$$
所以可以考虑对颜色进行分类，把颜色相同的物品放一起

1）对颜色分类，这样每一个颜色就是一个一个的段，先考虑`color[j-1] --> color[j]`，两个大段之间的 dp
然后再考虑段内元素之间的 dp

2）对于段内元素的 dp，记$$dp(i, 0/1, x)$$表示当前考虑段内的第$$i$$个位置，
`0/1`表示有没有选过当前段表示的颜色，`x`表示当前剩下的钱数

3）重点考虑段和段交接的时候，dp 的转移，假设当前段是$$j$$
如果`j = 0`，那么很显然`dp(0, 1, X-p[0]) = u[0] + K`（其中`p[0]，u[0]`分别表示第`0`个物品的价格和价值）

如果`j > 0`，那么`段首的状态 dp(0, ...)`，应该从上一个段末的状态转移过来，上一个段假设有`m+1`个元素（最后一个元素下标是`m`）
`for x = 枚举上一个状态还剩下多少钱`，有转移$$dp(0, 0, x) \leftarrow \max(last(m, 0, x), last(m, 1, x))$$
同时，$$dp(0, 1, x-p_0) \leftarrow \max(last(m, 0, x), last(m, 1, x)) + u_0 + K$$

4）这样每一段最开始（或者说段和段的交接）的状态我们都考虑到了，而段内的每个元素转移都是简单的线性 dp
`for i = 段内的每一个元素， for x = 枚举当前有多少钱`
转移不难，如果不选元素，$$dp(i, 0, x) \leftarrow dp(i-1, 0, x), dp(i, 1, x) \leftarrow dp(i-1, 1, x)$$
如果选元素，那么$$dp(i, 1, x) \leftarrow dp(i-1, 1, x+p_i) + u_i, \quad dp(i, 1, x) \leftarrow dp(i-1, 0, x+p_i) + u_i + K$$

考虑完最后一个段之后，对最后的段每个`dp(i, 0/1, x)`取一个`max`就是答案了

### G
**[G - Bar Cover](https://atcoder.jp/contests/abc383/tasks/abc383_g)**
你有一行$$n$$个单元格，第$$i$$个单元格的整数是$$A_i$$，另外你还有$$\lfloor \dfrac{N}{m} \rfloor$$块方格

每块方格覆盖了连续$$m$$个单元格，针对每一个$$i = 1, 2, \cdots, \lfloor \dfrac{N}{m} \rfloor$$段
求解，当恰好放置$$i$$个不重叠的单元的时候，被覆盖的部分的数字最大的和

#### 算法分析

> 注意到，$$m$$的上界只有$$5$$，能不能状态压缩？怎么状态压缩？但似乎不能做
这种一段一段的，可以考虑下区间 dp？那么区间 dp，它不是连续的段，状态怎么设计呢？
实际上，把空出的段也要设计到状态里面

记$$dp(l, r, x, y)$$，表示$$[l, r]$$区间，左边空出了$$x$$长度，右边空出了$$y$$长度
其中$$x \in [0, m), y \in [0, m)$$，（不能等于$$m$$，等于$$m$$的情况我们应该考虑**覆盖**）
然后考虑$$mid = (l+r)/2$$，区间的合并

1）中间**没有空位**，那么合并是简单的，$$dp(l,r,x,y) \leftarrow dp(l, mid, x, 0) + dp(mid+1, r, 0, y)$$
2）中间**需要覆盖**，中间至少需要$$\geqslant m$$的长度，假设左半边占了$$k$$的长度，右半边占了$$m - k$$的长度
（注意$$k \in [1, m)$$，两边的**空位**都不能是$$0$$， 为什么呢？）

>前面我们状态定义的时候，给两边预留的空位$$x, y \in [0, m)$$，最长的空位就是$$m-1$$了
那么左半部分，右端点如果空位是`0`，右半部分呢？左端点的空位就势必是`m`，这不符合状态定义

3）对于左半段$$[l,mid]$$来说，$$[l, l+x-1]$$这一段是空着的，从$$l+x$$开始**放东西**
那么左半段，最右边需要放$$k$$的长度，实际上**新放入的段**在$$[mid-k+1, mid]$$
所以有限制$$mid-k+1 \geqslant l+x$$

对于右半段呢？右半段保证$$[r-y+1, r]$$都是空着的，从$$r-y$$往前开始放东西
新放入的段的长度是$$m-k$$，$$[mid+1, mid+m-k]$$是新放入的段
那么有限制$$mid+m-k \leqslant r-y$$

4）考虑合并，`for k = [1, m-1]`，先必须满足上述限制，然后有转移
```math
dp(l, r, x, y) \leftarrow dp(l,mid,x, k) + dp(mid+1, r, m-x, y) + Mid
```
其中`Mid = 新放入的段的和`，新放入的段是$$[mid-k+1,  mid +m -k]$$，从而`Mid = S[mid+m-k] - S[mid-k]`

这样似乎就可以解决这个问题了

> 但是暴力的复杂度太高了，考虑归并，怎么归并呢？容易发现最优解具有凸性

记$$f(i)$$表示放$$i$$个数能够取到的最大的和，$$f(i+1), f(i+2)$$同理
注意到，放入$$i$$个数之后，最大的和能够在剩下的数中随便取，没有限制
而放入第$$i+1$$个数之后，就有限制了，最大的数被取走了，最好的情况只能取**次大数**
所以我们有$$f(i) = a, f(i+1) = a+b, f(i+2) = a+b+c, \quad a \geqslant b \geqslant c$$

而$$f(i) + f(i+2) = 2a+b+c$$，$$2f(i+1) = 2a + 2b$$，而$$b+c \leqslant 2b$$，实际上
我们有 $$f(i) + f(i+2) \leqslant 2f(i+1)$$，从而最优解有**上凸性**，换句话说，**差分数组是单调递减的**

实现的时候，我们用一个**凸包 Hull** 维护左半段和右半段
每一个 Hull 维护的是一个`vector`，`vector`存的实际上是**差分数组**，比如说，凸包有`2`个数
分别是第`1`段的最优解`f(1)`，第`2`段的最优解`f(2)`，那么`vector = { f(1), f(2)-f(1) }`
凸包的合并，就是对差分数组归并，并且按照**单调递减**顺序排列

用**闵可夫斯基和优化 dp**，就是在考虑`Hull(u) = max(Hull(u), Hull(v))`的时候
先根据差分数组，求前缀和算出每个`i`的实际值，然后对于每个`i`，`res(i) = max( U(i), V(i) )`
最后把`res`再转化成**差分数组**，返回差分数组即可

#### 算法实现
区间所有的线段，用凸包维护起来，凸包`Hull(u)`维护了区间`[l, r]`所有可能出现的状态
比如说，两侧分别有`3, 5`个空位，那么就对应凸包状态`Hull(u, 3, 5)`
```bash
using arr = array<int, 3>;

void solve() {
    int n, m;
    cin >> n >> m;
    vector<i64> a(n+1, 0);
    for (int i = 1; i <= n; i++) cin >> a[i];
    vector<i64> s(n+1, 0);
    for (int i = 1; i <= n; i++) s[i] = s[i-1] + a[i];

map<arr, vector<i64> > hull;
    int tot = 0;

auto _merge = [&](const vector<i64> &a, const vector<i64> &b) -> vector<i64> {
        vector<i64> res;
        if (a.size() == 0) return b;
        if (b.size() == 0) return a;

res.resize(a.size() + b.size());
        std::merge(a.begin(), a.end(), b.begin(), b.end(), res.begin(), greater<i64>());

return res;
    };

auto chkmax = [&](vector<i64> a, vector<i64> b) -> vector<i64> {
        if (a.size() == 0) return b;
        if (b.size() == 0) return a;
        vector<i64> res(max(a.size(), b.size()), -inf<i64>);

for (int i = 1; i < a.size(); i++) a[i] += a[i-1];
        for (int i = 1; i < b.size(); i++) b[i] += b[i-1];

for (int i = 0; i < res.size(); i++) {
            if (i < a.size()) chmax(res[i], a[i]);
            if (i < b.size()) chmax(res[i], b[i]);
        }
        for (int i = res.size()-1; i >= 1; i--) res[i] -= res[i-1];
        return res;
    };

auto solve = [&](auto &solve, int l, int r) -> int {
        int u = ++tot;
        if (l >= r) {
            if (m <= 1) hull[{u, 0, 0}].push_back(a[l]);
            return u;
        }
        int mid = (l + r) >> 1;
        int le = solve(solve, l, mid), ri = solve(solve, mid+1, r);

for (int x = 0; x < m; x++) for (int y = 0; y < m; y++) {
            if (x + y > r - l + 1) continue;
            hull[{u, x, y}] = _merge(hull[{le, x, 0}], hull[{ri, 0, y}]);

for (int k = 1; k < m; k++) {
                if (! (mid-k+1 >= l+x && mid+m-k <= r-y) ) continue;
                vector<i64> middle {s[mid+m-k] - s[mid-k]};
                auto side = _merge(hull[{le, x, k}], hull[{ri, m-k, y}]);
                auto now = _merge(side, middle);

hull[{u, x, y}] = chkmax(hull[{u, x, y}], now);
            }
        }

return u;
    };

auto id = solve(solve, 1, n);
    i64 sum = 0;
    for (auto x : hull[{id, 0, 0}]) {
        sum += x;
        cout << sum << " ";
    }
    cout << "\n";
}
```

## 闵可夫斯基和
### 凸集和闵可夫斯基和
**凸集的定义**
对于定义在凸集$$\Omega \subset \mathbb{R}^n$$上的函数$$f:\Omega \to \mathbb{R}$$，$$f$$是凸函数当且仅当
对于任意$$\boldsymbol{x, y} \in \Omega$$和任意的$$\alpha \in (0, 1)$$，都有
```math
\displaystyle
f(\alpha \boldsymbol{x} + (1-\alpha) \boldsymbol{y}) \leqslant \alpha f(\boldsymbol{x}) + (1-\alpha)f(\boldsymbol{y})
```
实际上这是关于下凸的定义，如果是上凸，那么有
```math
\displaystyle
\alpha f(\boldsymbol{x}) + (1-\alpha)f(\boldsymbol{y}) \leqslant f(\alpha \boldsymbol{x} + (1-\alpha) \boldsymbol{y})
```
一般情况，在算法竞赛问题中，尝尝取$$\alpha = \dfrac{1}{2}$$

**闵可夫斯基和**
点集$$P$$和点集$$Q$$的闵可夫斯基和定义为$$P+Q = \\{a+b | a\in P, b \in Q\\}$$
把点集$$Q$$中的每一个点，看作一个向量（比如点$$E$$，那么就是$$\overrightarrow{OE}$$）
将点集$$P$$中的每一个点，沿着这些向量**平移**，最终得到的结果的集合，就是点集$$P+Q$$

算法实现过程如下图所示
![](/media/content_img/minkowski01.png)
![](/media/content_img/minkowski02.png)

**凸集的一些性质**
> 若点集$$P, Q$$是凸集，那么闵可夫斯基和$$P+Q$$也是凸集

证明，不妨设$$e, f \in P+Q$$，同时我们有$$a, b \in P, c, d\in Q$$，并且$$e = a+c, f = b + d$$，对于$$\forall t \in (0, 1)$$，有
```math
te+(1-t)f = t(a+c)+(1-t)(b+d) 
\\= (ta+(1-t)b) + (tc + (1-t)d) \in P + Q
```

> 若点集$$P, Q$$为凸集，则其闵可夫斯基和$$P+Q$$的边集，是由凸集$$P, Q$$的边，按极角排序后连接的结果

证明，如图所示，可以旋转坐标轴，令$$\boldsymbol{x}$$和$$P$$的一边重合，一开始的时候
我们选择**重合$$P$$最下方的那条边**
接着执行向量运算的时候，我们找到$$Q$$最低的点（图中所示）$$F$$，执行向量和$$\boldsymbol{x} + \overrightarrow{OF} = \overrightarrow{OC}$$
不难发现，$$C$$一定位于凸包$$P+Q$$的边界上
另外$$\overrightarrow{CC'}$$可以看作是**向量平行四边形的边**$$\overrightarrow{OF}$$**平移**得到
这样我们就证明了，以**$$P$$最低点的边**，和$$Q$$最低点的边，恰好构成了**凸包第一段转折**
即$$\overline{AC}, \overline{CC'}$$

接着考虑扔掉上述**构造好的两条边**，继续让$$\boldsymbol{x}$$重合$$P$$的最低点
根据数学归纳法，同理可以得到**第二段转折**

实际上，$$P+Q$$的边，按**从左到右，极角序，顺时针顺位排下来**

### 算法实现
将凸集极角排序（实际上，叉积，顺序$$a < b$$，当且仅当$$\boldsymbol{a} \times \boldsymbol{b} < 0$$）
将$$P_0 + Q_0$$作为并集的第一个元素，然后对`vector`归并，依次放边即可
```bash
template <typename T>
struct Point {
    T x, y;
    Point(T x = 0, T y = 0) : x(x), y(y) {}

friend Point operator+ (const Point &a, const Point &b) {
        return {a.x + b.x, a.y + b.y};
    }
    friend Point operator- (const Point &a, const Point &b) {
        return {a.x - b.x, a.y - b.y};
    }
    // dot
    friend T operator* (const Point &a, const Point &b) {
        return a.x * b.x + a.y * b.y;
    }
    // cross
    friend T operator^ (const Point &a, const Point &b) {
        return a.x * b.y - a.y * b.x;
    }
};

template <typename T>
vector<Point<T> > minkowski(vector<Point<T> > a, vector<Point<T> > b) {
    vector<Point<T> > res{a[0] + b[0]};
    // diffrence
    for (int i = 0; i + 1 < (int)a.size(); i++) a[i] = a[i+1] - a[i];
    for (int i = 0; i + 1 < (int)b.size(); i++) b[i] = b[i+1] - b[i];
    a.pop_back(), b.pop_back();

res.resize(a.size() + b.size() + 1);
    merge(a.begin(), a.end(), b.begin(), b.end(), c.begin()+1, [](const Point<T> &a, const Point<T> &b) {
        return (a ^ b) < 0;
    });

for (int i = 1; i < res.size(); i++) res[i] = res[i] + res[i-1];
    return res;
}
```

## 凸集的保凸操作
**函数的加权和**
![](/media/content_img/convex01.png)

**函数的取 max**
![](/media/content_img/convex02.png)

**函数的合成，凸集的交**
![](/media/content_img/convex03.png)

## 扩展欧几里得算法
### exgcd
注意到 exgcd 的原理并不难，实际上对于$$a, b$$，我们可以求解$$x, y$$
满足$$ax + by = \text{gcd}(a, b)$$
```math
\begin{cases}
ax + by = \text{gcd}(a, b) \\
by' + (a \bmod b)x' = \text{gcd}(b, a \bmod b)
\end{cases}
```
从而有递归$$\text{exgcd}(a, b, x, y) \to \text{exgcd}(b, a\bmod b, y', x')$$
而$$by'+ (a - \lfloor a / b \rfloor \cdot b) x' = \text{gcd}(b, a\bmod b)  = \text{gcd}(a, b)$$

写成关于$$a, b$$的，$$(y' - \lfloor a/b \rfloor x') \cdot b + ax' = \text{gcd}(a, b)$$
那么我们递归求解$$x', y'$$之后，令$$x \leftarrow x', \quad y \leftarrow y' - \lfloor a/b \rfloor \cdot x'$$

> $$(a, b) = 1 \Leftrightarrow \exists x, y, ax+by = 1$$

### 求一般线性丢番图方程
根据扩展欧几里得算法，求解一般线性丢番图方程$$ax+by = c$$
1）先用扩展欧几里得算法求出$$g = \text{gcd}(a, b)$$，如果$$g \nmid c$$，无解
2）扩展欧几里得可以求出**特解**$$x_0, y_0$$，即满足$$ax_0 + by_0 = g$$

那么关于$$ax + by = c$$的解是多少呢？很明显，我们有$$x = x_0 \cdot c/g, \ y = y_0 \cdot c/g$$

实际上我们可以写成$$a\left( \dfrac{x_0}{g} \right) c + b \left( \dfrac{y_0}{g} \right) c = c$$

3）对于一般线性丢番图方程$$ax + by = c$$，用 exgcd 求出特解$$x = x_0 \cdot \dfrac{c}{g}, \ y = y_0 \cdot \dfrac{c}{g}$$，不难发现
```math
\displaystyle
a\left(x + k \dfrac{b}{\text{gcd}(a, b)} \right) + b\left(x - k \dfrac{a}{\text{gcd}(a, b)} \right) = c
```
对于$$k \in \mathbb{Z}$$，仍然满足，通解是
```math
\displaystyle
\begin{cases}
x + k \dfrac{b}{\text{gcd}(a, b)} \\
y - k \dfrac{a}{\text{gcd}(a, b)}
\end{cases}
```

> 实际应用中，经常要求关于$$ax+by = c$$的最小正整数解

1）先求$$ax+by = g = \text{gcd}(a, b)$$，exgcd 可以同时求出一组特解$$x_0, y_0$$

2）令$$x_0 = x_0 \cdot \dfrac{c}{g}, \quad y_0 = y_0 \cdot \dfrac{c}{g}$$

3）考虑$$ax+by = c$$的通解，$$x = x_0 +\dfrac{b}{g}k$$，$$minx = x_0 \bmod \dfrac{b}{g}$$
（这里用 safeMod）

**[青蛙的约会](https://ac.nowcoder.com/acm/contest/21094/A)**
```bash
constexpr i64 safeMod(i64 x, i64 m) {
    x %= m;
    if (x < 0) {
        x += m;
    }
    return x;
}

constexpr std::array<i64, 3> exgcd(i64 a, i64 b) {
    if (b == 0) {
        return {a, 1, 0};
    }
    auto [g, x, y] = exgcd(b, a % b);
    return {g, y, x - a/b * y};
}

void solve() {
    i64 x, y, m, n, L;
    cin >> x >> y >> m >> n >> L;

i64 a = m - n, c = y - x;
    if (a < 0) a = -a, c = -c;

auto [g, u, v] = exgcd(a, L);
    if (c % g) cout << "Impossible\n";
    else {
        u *= c / g, v *= c / g;
        u = safeMod(u, L/g);
        cout << u << "\n";
    }
}
```

## 模数和逆元
### exgcd 求逆元
求$$a \equiv a^{-1} (\bmod p)$$，这里不要求$$p$$是素数，但是必须$$(a, p) = 1$$
我们有$$ax + py = 1$$，等价于$$ax \equiv 1 (\bmod p)$$，即$$x \equiv a^{-1} (\bmod p)$$

用 exgcd 算法求出满足条件的**最小正整数**$$x$$即可
```bash
auto [g, x] = exgcd(a, mod, x, y);
if (x < 0) x = safeMod(x, mod);
```

### Barrett 模乘
要求$$a \cdot b \bmod m$$，也就是我们要找到余数$$v$$，$$0 \leqslant v < m$$，其中$$a /b$$商是$$x$$
```math
v = a \cdot b - x\cdot m \quad (0 \leqslant v < m)
```
另外，$$x = \lfloor \dfrac{a \cdot b}{m} \rfloor = \lfloor \dfrac{a \cdot b}{2^k} \cdot \dfrac{2^k}{m} \rfloor$$，这里$$k$$可以取很大的数，比如`unsigned long long 的 64`

记$$im = \lfloor \dfrac{2^k}{m} \rfloor$$，取$$x$$的近似，$$x' = \lfloor \dfrac{a \cdot b}{2^k} \cdot im \rfloor$$

这里取`k = 64 (对应 unsigned long long 的位数)`，那么$$im = \lfloor \dfrac{2^k}{m} \rfloor = \dfrac{2^k + m - 1}{m}= \dfrac{2^{64} - 1}{m} + 1$$
所以可以构造`im = (u64)(-1) / m + 1`

另外，$$x^{\prime} = \lfloor \dfrac{a\cdot b \cdot im}{2^k} \rfloor$$，实际上
`u64 x = u64( (u128)(a) * (u128)(b) * im ) >> 64`，这样就求解出了$$x^{\prime}$$

> 有一个结论，近似的$$x^{\prime}$$和$$x$$之间，最多相差$$1$$

也就是说，$$x^{\prime} = x$$，或者$$x^{\prime} = x - 1$$，那么最终的$$v$$呢
我们有$$v = a \cdot b - x\cdot m$$，或者$$v = a \cdot b - (x-1)\cdot m = a\cdot b - x\cdot m + m$$

换句话说，我们根据$$v = a\cdot b - x \cdot m$$求出$$v$$，如果发现$$v \geqslant m$$
直接$$v += m$$即可

（注意，这里的$$v$$用的是**无符号整型**，如果发现$$v \geqslant m$$，说明$$v$$太大了，负数，溢出了）
（所以这里令`v += m`）

```bash
template<u32 P>
using ModInt = ModIntBase<u32, P>;
template<u64 P>
using ModInt64 = ModIntBase<u64, P>;

struct Barrett {
public:
    Barrett(u32 m_) : m(m_), im((u64)(-1) / m_ + 1) {}

constexpr u32 mod() const {
        return m;
    }

constexpr u32 mul(u32 a, u32 b) const {
        u64 z = a;
        z *= b;
        
        u64 x = u64((u128(z) * im) >> 64);
        
        u32 v = u32(z - x * m);
        if (m <= v) {
            v += m;
        }
        return v;
    }

private:
    u32 m;
    u64 im;
};
```

### 逆元表
可以根据线性递推来求，一开始的时候，我们有$$1^{-1} \equiv 1(\bmod p)$$，现在要打表求$$i^{-1}(\bmod p)$$

1）不妨设$$p = ki + j$$，其中$$k = \lfloor p / i \rfloor$$，$$j = p \bmod i, \quad 0 \leqslant j < i$$
2）同时注意到，条件等价于$$ki + j \equiv 0(\bmod p)$$，两边同时$$\times i^{-1}j^{-1}$$
我们有$$i^{-1} \equiv -kj^{-1} (\bmod p)$$
3）将$$k, j$$代入，我们有
```math
i^{-1} \equiv \begin{cases} 1 & (i = 1) \\ -\lfloor \dfrac{p}{i} \rfloor (p \bmod i)^{-1} \end{cases}
```
递推就是
```bash
inv[i] = (P - P/i) * inv[P % i]
```

```bash
using Z = ModInt<998244353>;

template<u32 P>
struct InvTable {
    int _n;
    std::vector<Z> _inv;

InvTable() : _n(0), _inv{1} {};
    InvTable(int n) : InvTable() {
        init(n);
    }

void init(int n) {
        if (n <= _n) return;
        _inv.resize(n+1);

_inv[1] = 1;
        for (int i = 2; i <= n; i++) {
            _inv[i] = Z(P - P / i) * _inv[P % i];
        }

_n = n;
    }
};
```

### 阶乘递推
求组合数的时候，经常涉及到阶乘的逆元，阶乘的递推，其实并不复杂
1）首先，我们有$$x! = (x-1)! \cdot x$$
```bash
_fac[i] = _fac[i-1] * i
```

2）关于阶乘的逆元，也有递推，$$\dfrac{1}{(m-1)!} = \dfrac{1}{m!} \cdot m$$
顺便还可以求出逆元表$$\dfrac{1}{m} = \dfrac{1}{m!} \cdot (m-1)!$$
```bash
_invfac[m] = _fac[m].inv();
for (int i = m; i > n; i--) {
	_invfac[i - 1] = _invfac[i] * i;
	_inv[i] = _invfac[i] * _fac[i - 1];
}
```

上一篇：2024年12月一些算法比赛（一）下一篇：图论遗留问题，最短路，二分图

看完文章有啥想法