BITACM 2025 暑期集训 B组笔记

Notes

Notes

[2025-07-15] DP基础

0/1背包

有$n\in\mathbb N$个物品，重量分别为$\vec w\in(\N^*)^n$，价值分别为$\vec v\in(\N^*)^n$。背包容量为$c$，求包里能装下物品的最大价值。

设$f_{i,j}$表示对于前$i$个物品，容量为$j$的背包能装下物品的最大价值。

递推边界：

\[f_{0,i}=0\ (i\in\N,\ i\leq c) \]

状态转移方程：

\[f_{i+1,j}= \left\{ \begin{align*} &0&&(j<w_i)\\ &\max\left(f_{i,j},f_{i,j-w_i}+v_i\right)&&(j\geq w_i) \end{align*} \right. \ (i,j\in\N,\ i<n,\ j\leq c) \]

时间复杂度：$O(n\cdot c)$。
空间复杂度：$O(n\cdot c)$，采用滚动数组优化后可得$O(c)$。

完全背包

有$n\in\mathbb N$种物品，重量分别为$\vec w\in(\N^*)^n$，价值分别为$\vec v\in(\N^*)^n$，每种物品可以取任意个。背包容量为$c$，求包里能装下物品的最大价值。

设$f_i$表示容量为$i$的背包能装下物品的最大价值。

状态转移方程：

\[f_i= \left\{ \begin{align*} &0&&\left(i<\min\vec w\right)\\ &\max_{j\in\N,\ j<n,\ w_j\leq i}\left(f_{i-w_j}+v_j\right)&&\left(i\geq\min\vec w\right) \end{align*} \right. \ (i\in\N,\ i\leq c) \]

时间复杂度：$O(c\cdot n)$。
空间复杂度：$O(c)$。

多重背包

有$n\in\mathbb N$种物品，重量分别为$\vec w\in(\N^*)^n$，价值分别为$\vec v\in(\N^*)^n$，最多分别可以取$\vec m\in(\N^*)^n$个。背包容量为$c$，求包里能装下物品的最大价值。

设$f_{i,j}$表示对于前$i$种物品，容量为$j$的背包能装下物品的最大价值。

递推边界：

\[f_{0,i}=0\ (i\in\N,\ i\leq c) \]

状态转移方程：

\[f_{i+1,j}= \left\{ \begin{align*} &0&&(j<w_i)\\ &\max_{k=0}^{m_i}\left(f_{i,j-k\cdot w_i}+k\cdot v_i\right)&&(j\geq w_i) \end{align*} \right. \ (i,j\in\N,\ i<n,\ j\leq c) \]

时间复杂度：$O(c\sum\vec m)$。
空间复杂度：$O(n\cdot c)$，采用滚动数组优化后可得$O(c)$。

优化

对于第$i$种物品，分解$m_i=\sum_{j=0}^{p_i-1}2^{j}+(m_i-(2^{p_i}-1))$，其中$p_i=\left\lfloor\log(m_i+1)\right\rfloor$。分组使第$j\ (j\in\mathbb N,\ j<p_i)$组有$2^j$个物品，第$p_i$组有$m_i-(2^{p_i}-1)$个物品。

对分组后的物品应用上述朴素方法。

时间复杂度：$O(c\sum_{i=0}^{n-1}\log m_i)$。
空间复杂度：$O(n\cdot c+\sum_{i=0}^{n-1}\log m_i)$，采用滚动数组优化后可得$O(c+\sum_{i=0}^{n-1}\log m_i)$。

最长上升子序列（LIS）

用$min\_end_i$维护长为$(i+1)$的末尾元素最小的上升子序列的末尾元素位置，因为一个上升子序列的末尾元素越小，之后就越容易延长。如果当前元素可以取代$min\_end$中某一项而不破坏其顺序（即对于当前的$a_i$，$(\exist j\in\mathbb N,\ j<m)\ ((j=0\vee a_{min\_end_{j-1}}<a_i)\wedge(j=m-1\vee a_i<a_{min\_end_{j+1}}))$，其中$m$是$min\_end$的长度），则用当前元素替换掉那一项；如果当前元素大于min_end中所有元素，则将当前元素插入min_end末尾。这样，min_end.size()就始终是当前能得到的LIS长度。

不过，此时$min\_end$并非我们想要的，最终LIS各元素的位置，因为可能在我们得到了LIS后，某些元素又被后来的更小的元素覆盖了。不过，$min\_end$最后那个元素一定对应LIS的最后一个元素。于是，我们可以记录每个元素放入$min\_end$时，它的前一个元素是哪个；从$min\_end_{m-1}$开始沿着前驱走，就可以复原LIS的各个元素了。

/**
 * @return indices of elements consisting
 * 	one longest increasing subsequence of @param arr
 */
template <typename T, typename Cmp = std::less<T>>
std::vector<size_t> longestIncrSubseq(const std::vector<T> &arr, Cmp cmp = Cmp()) {
	auto cmpVal = [&](size_t lhs, size_t rhs) { return cmp(arr[lhs], arr[rhs]); };
	std::vector<size_t> min_end, res;
	auto pre = std::vector<size_t>(arr.size());
	for (size_t i = 0; i != arr.size(); ++i) {
		size_t j = std::lower_bound(min_end.begin(), min_end.end(), i, cmpVal) -
				   min_end.begin();
		if (j == min_end.size()) {
			min_end.push_back(i);
		} else {
			min_end[j] = i;
		}
		pre[i] = (j ? min_end[j - 1] : size_t(-1));
	}
	res.reserve(min_end.size());
	for (size_t p = min_end.back(); ~p; p = pre[p]) {
		res.push_back(p);
	}
	std::reverse(res.begin(), res.end());
	return res;
}

[2025-07-16] 树状数组和线段树

树状数组

用于在单点修改时维护一个数组的区间和（通常只需维护前缀和）。

懒得抽象了。这里就只写最简单的，维护数的加法的树状数组罢。理论上只需要一个运算满足结合律即可用树状数组维护，但幺元、交换律会让实现更加简单，有时也可以让效率更高。

设原数组为$\vec a\in\mathbb F^n\ (n\in\mathbb N)$。树状数组维护一个$\vec t\in\mathbb F^{n+1}$使得

\[t_i=\sum_{j=i-\mathrm{lsb}(i)}^{i-1}a_i\ (i\in\N,\ i\leq n) \]

其中$\mathrm{lsb}$是指Lowest Significant Bit，更常见的名字是$\mathrm{lowbit}$。

详见OI-Wiki。

下面是一个功能比较完备的树状数组模板，但它对于数以外的对象和运算（例如矩阵乘法、字符串的拼接等）支持很差。将幺元作为模板参数并不是合适的做法。

template <typename Iter>
using RequireInputIter =
	std::enable_if_t<
		std::is_convertible<
			typename std::iterator_traits<Iter>::iterator_category,
			std::input_iterator_tag>::value>;

template <typename T, typename Oper = std::plus<T>, T id_elem = T()>
class FenwickTree {
public:
	inline static size_t lowbit(size_t x) { return (x & (-x)); }

	inline explicit FenwickTree(size_t n = 0) : m_tree(n + 1, id_elem) {}
	inline explicit FenwickTree(std::initializer_list<T> list) : m_tree(list) {
		m_build();
	}
	inline FenwickTree(size_t n, const T &value) : m_tree(n + 1, value) {
		m_tree.front() = id_elem;
		m_build();
	}
	template <typename Iter, typename = RequireInputIter<Iter>>
	inline FenwickTree(Iter begin, Iter end)
		: m_tree(std::distance(begin, end) + 1) {
		m_tree.front() = id_elem;
		std::copy(begin, end, m_tree.begin() + 1);
		m_build();
	}

	inline size_t treeSize() const { return m_tree.size(); }
	inline size_t size() const { return (m_tree.size() - 1); }

	inline void resize(size_t n, const T &val = id_elem) {
		if ((++n) <= m_tree.size()) {
			m_tree.resize(n);
		} else {
			size_t old_sz = m_tree.size();
			m_tree.resize(n, val);
			m_rebuild(old_sz);
		}
	}

	inline void assign(std::initializer_list<T> list) {
		assign(list.begin(), list.end());
	}
	inline void assign(size_t n = 0) { m_tree.assign(n + 1, id_elem); }
	inline void assign(size_t n, const T &value) {
		m_tree.assign(n + 1, value);
		m_tree.front() = id_elem;
		m_build();
	}
	template <typename Iter, typename = RequireInputIter<Iter>>
	inline void assign(Iter begin, Iter end) {
		m_tree.resize(std::distance(begin, end) + 1);
		std::copy(begin, end, m_tree.begin() + 1);
		m_build();
	}

	/**
	 * @brief add diff to the element at index
	 */
	inline void modify(size_t index, const T &diff) {
		m_range_check(index);
		for (++index; index < m_tree.size(); index += lowbit(index)) {
			m_tree[index] = m_oper(m_tree[index], diff);
		}
	}

	/**
	 * @return the sum of [0, min(n, size()))
	 */
	inline T query(size_t n = -1) const {
		if (n >= m_tree.size()) {
			n = m_tree.size() - 1;
		}
		T res = id_elem;
		for (; n; n -= lowbit(n)) {
			res = m_oper(res, m_tree[n]);
		}
		return res;
	}
	inline T operator[](size_t n) const { return query(n); }

protected:
	std::vector<T> m_tree; // m_tree[i] maintains the sum of data[i - lowbit(i), i)
	Oper m_oper;

	inline void m_range_check(size_t index) const {
		if (index + 1 >= m_tree.size()) {
			std::__throw_out_of_range_fmt(__N("FenwickTree::__range_check: index "
											  "(which is %zu) >= this->size() "
											  "(which is %zu)"),
										  index, m_tree.size() - 1);
		}
	}
	inline void m_build() {
		for (size_t i = 1, j; i < m_tree.size(); ++i) {
			j = i + lowbit(i);
			if (j < m_tree.size()) {
				m_tree[j] = m_oper(m_tree[j], m_tree[i]);
			}
		}
	}
	inline void m_rebuild(size_t old_tree_sz) {
		for (size_t i = 1, j; i < m_tree.size(); ++i) {
			j = i + lowbit(i);
			if (j >= old_tree_sz && j < m_tree.size()) {
				m_tree[j] = m_oper(m_tree[j], m_tree[i]);
			}
		}
	}
};

时间复杂度：修改和查询都是$O(\log n)$。
空间复杂度：$O(n)$。

线段树

[2025-07-17] DFS、BFS、最短路和拓扑排序

贴几个模板罢。

最短路

template <typename Weight>
inline std::vector<std::vector<Weight>>
floyd(std::vector<std::vector<Weight>> weights) {
	size_t n = weights.size();
	if (!n) {
		return {};
	}
	assert(n == weights[0].size());
	for (size_t i = 0, j, k; i != n; ++i) {
		for (j = 0; j != n; ++j) {
			for (k = 0; k != n; ++k) {
				if (weights[j][i] != std::numeric_limits<Weight>::max() &&
					weights[i][k] != std::numeric_limits<Weight>::max()) {
					weights[j][k] = std::min(weights[j][k], weights[j][i] + weights[i][k]);
				}
			}
		}
	}
	return weights;
}

template <typename Weight>
inline std::vector<Weight>
dijkstra(const std::vector<std::vector<std::pair<size_t, Weight>>> &adj, // {to, weight}
		 size_t src) {
	size_t n = adj.size();
	if (!n) {
		return {};
	}
	using Adj = std::pair<Weight, size_t>; // {weight, to}
	std::priority_queue<Adj, std::vector<Adj>, std::greater<Adj>> pq;
	auto dist = std::vector<Weight>(n, std::numeric_limits<Weight>::max());
	pq.emplace(0, src);
	dist[src] = 0;
	auto vis = std::vector<bool>(n);
	while (pq.size()) {
		auto u = pq.top().second;
		pq.pop();
		if (vis[u]) {
			continue;
		}
		vis[u] = true;
		for (auto &pr : adj[u]) {
			auto &v = pr.first;
			auto &w = pr.second;
			if (dist[v] > dist[u] + w) {
				pq.emplace(dist[v] = dist[u] + w, v);
			}
		}
	}
	return dist;
}

template <typename Weight>
inline std::vector<Weight>
bellmanFord(size_t n,
			const std::vector<std::tuple<size_t, size_t, Weight>> &edges, // {from, to, weight}
			size_t start) {
	if (!n) {
		return {};
	}
	auto dist = std::vector<Weight>(n, std::numeric_limits<Weight>::max());
	dist[start] = 0;
	for (size_t cnt = 0; cnt != n; ++cnt) {
		bool flag = true;
		for (auto &edge : edges) {
			auto &u = std::get<0>(edge), &v = std::get<1>(edge);
			auto &w = std::get<2>(edge);
			if (dist[u] != std::numeric_limits<Weight>::max() &&
				dist[v] > dist[u] + w) {
				dist[v] = dist[u] + w;
				flag = false;
			}
		}
		if (flag) {
			return dist;
		}
	}
	return {}; // A circle of neg-weight edges can be reached from the start point.
}

template <typename Weight>
inline std::vector<Weight>
spfa(const std::vector<std::vector<std::pair<size_t, Weight>>> &adj, // {to, weight}
	 size_t start) {
	size_t n = adj.size();
	if (!n) {
		return {};
	}
	auto dist = std::vector<Weight>(n, std::numeric_limits<Weight>::max());
	auto inq = std::vector<bool>(n);
	auto cnt = std::vector<size_t>(n);
	dist[start] = 0;
	std::queue<size_t> q;
	q.emplace(start);
	while (q.size()) {
		auto u = q.front();
		inq[u] = false;
		q.pop();
		for (auto &e : adj[u]) {
			auto v = e.first;
			auto w = e.second;
			if (dist[v] > dist[u] + w) {
				dist[v] = dist[u] + w;
				/**
				 * the shortest path between 2 vertices consists of at most
				 * (n - 1) edgess
				 */
				if ((cnt[v] = cnt[u] + 1) >= n) {
					return {};
				}
				if (!inq[v]) {
					inq[v] = true;
					q.emplace(v);
				}
			}
		}
	}
	return dist;
}

拓扑排序

inline std::vector<size_t>
toposort(const std::vector<std::vector<size_t>> &adj_vers) {
	size_t n = adj_vers.size();
	auto indeg = std::vector<uint32_t>(n);
	for (size_t from = 0; from != n; ++from) {
		for (size_t to : adj_vers[from]) {
			++indeg[to];
		}
	}
	std::vector<size_t> res, zero_indeg_vers;
	for (size_t i = 0; i != n; ++i) {
		if (!indeg[i]) {
			zero_indeg_vers.push_back(i);
		}
	}
	while (zero_indeg_vers.size()) {
		size_t from = zero_indeg_vers.back();
		zero_indeg_vers.pop_back();
		res.push_back(from);
		for (size_t to : adj_vers[from]) {
			if (!(--indeg[to])) {
				zero_indeg_vers.push_back(to);
			}
		}
	}
	return ((res.size() == n) ? res : std::vector<size_t>());
}

[2025-07-19] 单调栈和单调队列

单调栈

单调队列

template <typename T, typename Cmp = std::less<T>, typename Dq = std::deque<T>>
class MonoDeque {
public:
	inline MonoDeque(const Cmp &cmp = Cmp()) : m_cmp(cmp) {}

	inline size_t size() const noexcept { return m_dq.size(); }
	inline bool empty() const noexcept { return m_dq.empty(); }

	inline const T &front() const { return m_dq.front(); }
	inline const T &back() const { return m_dq.back(); }

	inline void pushFront(const T &val) {
		popFrontFor(val);
		m_dq.push_front(val);
	}
	inline void pushFront(T &&val) {
		popFrontFor(val);
		m_dq.push_front(std::move(val));
	}
	template <typename... Args>
	inline void emplaceFront(Args &&...args) {
		auto t = T(std::forward(args)...);
		popFrontFor(t);
		m_dq.push_front(std::move(t));
	}

	inline void pushBack(const T &val) {
		popBackFor(val);
		m_dq.push_back(val);
	}
	inline void pushBack(T &&val) {
		popBackFor(val);
		m_dq.push_back(std::move(val));
	}
	template <typename... Args>
	inline void emplaceBack(Args &&...args) {
		auto t = T(std::forward<Args>(args)...);
		popBackFor(t);
		m_dq.push_back(std::move(t));
	}

	inline void popFront() { m_dq.pop_front(); }
	inline void popFrontFor(const T &val) {
		while (m_dq.size() && !m_cmp(val, m_dq.front())) {
			popFront();
		}
	}

	inline void popBack() { m_dq.pop_back(); }
	inline void popBackFor(const T &val) {
		while (size() && !m_cmp(m_dq.back(), val)) {
			popBack();
		}
	}

	inline void clear() noexcept { m_dq.clear(); }

	inline void swap(MonoDeque &another) noexcept {
		std::swap(m_cmp, another.m_cmp);
		std::swap(m_dq, another.m_dq);
	}

protected:
	Cmp m_cmp;
	Dq m_dq;
};

[2025-07-20] 贪心、二分和双指针

贪心

二分

模板。

/**
 * @brief behaves as if it generates func[beg, end) and
 * 	performs std::lower_bound on it
 */
template <typename T, typename Ret, typename Func,
		  typename Cmp = std::less<Ret>>
inline T lowerBound(T beg, T end, const Ret &val, Func func = Func(), Cmp cmp = Cmp()) {
	while (beg < end) {
		T mid = beg + (end - beg) / 2;
		if (cmp(func(mid), val)) {
			beg = mid + 1;
		} else {
			end = mid;
		}
	}
	return beg;
}

/**
 * @brief behaves as if it generates func[beg, end) and
 * 	performs std::upper_bound on it
 */
template <typename T, typename Ret, typename Func,
		  typename Cmp = std::less<Ret>>
inline T upperBound(T beg, T end, const Ret &val, Func func = Func(), Cmp cmp = Cmp()) {
	while (beg < end) {
		T mid = beg + (end - beg) / 2;
		if (cmp(val, func(mid))) {
			end = mid;
		} else {
			beg = mid + 1;
		}
	}
	return beg;
}

双指针

[2025-07-23] 位运算和状压DP

位运算

状压DP

最短哈密顿回路

设$G=(V,E)$是一个图，若$G$中一条路径通过且仅通过每一个顶点一次，称这条路径为哈密顿路径（Hamiltonian path）。若$G$中一个回路通过且仅通过每一个顶点一次，称这个环为哈密顿回路（Hamiltonian cycle）。若一个图存在哈密顿回路，就称为哈密顿图（Hamiltonian/traceable graph）。

/**
 * @return the shortest Hamiltonian cycle beginning and ending at vertex 0
 */
template <typename Weight>
inline std::vector<size_t> shortestHamiltonianCycle(
	const std::vector<std::vector<Weight>> &adj_mat,
	Weight inf = std::numeric_limits<Weight>::max()) {
	size_t n = adj_mat.size();
	if (!n) {
		return {};
	}
	if (n == 1) {
		return {0};
	}
	auto dp = std::vector<std::vector<Weight>>(uint64_t(1) << n,
											   std::vector<Weight>(n, inf));
	auto pre = std::vector<std::vector<size_t>>(
		uint64_t(1) << n,
		std::vector<size_t>(n, size_t(-1)));
	dp[1][0] = 0;
	for (uint64_t set = 1; set < (uint64_t(1) << n); set += 2) {
		for (size_t cur = (set != 1); cur < n; ++cur) {
			if (!(set & (1 << cur))) {
				continue;
			}
			for (size_t nxt = 1; nxt < n; ++nxt) {
				if ((set & (1 << nxt)) || adj_mat[cur][nxt] == inf) {
					continue;
				}
				uint64_t nxt_set = set | (1 << nxt);
				if (dp[set][cur] + adj_mat[cur][nxt] < dp[nxt_set][nxt]) {
					pre[nxt_set][nxt] = cur;
					dp[nxt_set][nxt] = dp[set][cur] + adj_mat[cur][nxt];
				}
			}
		}
	}
	uint64_t full_set = (uint64_t(1) << n) - 1;
	size_t last = 1;
	for (size_t i = 2; i < n; ++i) {
		if (adj_mat[i][0] != inf &&
			dp[full_set][i] + adj_mat[i][0] <
				dp[full_set][last] + adj_mat[last][0]) {
			last = i;
		}
	}
	std::vector<size_t> res;
	res.reserve(n + 1);
	while (full_set) {
		res.emplace_back(last);
		size_t tmp = last;
		last = pre[full_set][last];
		full_set ^= (1 << tmp);
	}
	std::reverse(res.begin(), res.end());
	res.emplace_back(0);
	return res;
}

[2025-07-26] 稀疏表、左偏树和LCA

稀疏表

template <typename Iter>
using RequireInputIter =
	std::enable_if_t<
		std::is_convertible<
			typename std::iterator_traits<Iter>::iterator_category,
			std::input_iterator_tag>::value>;

template <typename T, typename Cmp = std::less<>>
struct Min {
	inline const T &operator()(const T &lhs, const T &rhs) const {
		return std::min(lhs, rhs, Cmp());
	}
};

template <typename T, typename Cmp = std::less<>>
struct Max {
	inline const T &operator()(const T &lhs, const T &rhs) const {
		return std::max(lhs, rhs, Cmp());
	}
};

template <typename _T, typename _Oper = Min<_T>>
class SparseTable {
public:
	using Oper = _Oper;
	using Elem = _T;
	using Table = std::vector<std::vector<Elem>>;

	inline static unsigned floorLogn2(uint64_t x) {
#if __cplusplus >= 202002L
		return (std::bit_width(x) - 1);
#else
		return (x ? (sizeof(uint64_t) * 8 - __builtin_clzll(x) - 1) : -1);
#endif
	}

	SparseTable(const Oper &oper = Oper()) : m_oper(oper) {}
	template <typename Iter, typename = RequireInputIter<Iter>>
	SparseTable(Iter begin, Iter end, const Oper &oper = Oper())
		: m_oper(oper), m_data(std::distance(begin, end)) {
		for (size_t i = 0; i != size(); ++i) {
			m_data[i].resize(floorLogn2(size() - i) + 1);
			m_data[i][0] = *(begin++);
		}
		for (size_t i = 1; (size_t(1) << i - 1) < size(); ++i) {
			for (size_t j = 0; j + (size_t(1) << i) <= size(); ++j) {
				m_data[j][i] = m_oper(m_data[j][i - 1],
									  m_data[j + (size_t(1) << (i - 1))][i - 1]);
			}
		}
	}

	inline void assign(const Oper &oper = Oper()) {
		m_oper = oper;
		m_data.clear();
	}
	template <typename Iter, typename = RequireInputIter<Iter>>
	inline void assign(Iter begin, Iter end, const Oper &oper = Oper()) {
		m_oper = oper;
		m_data.resize(std::distance(begin, end));
		for (size_t i = 0; i != size(); ++i) {
			m_data[i].resize(floorLogn2(size() - i) + 1);
			m_data[i][0] = *(begin++);
		}
		for (size_t i = 1; (size_t(1) << (i - 1)) < size(); ++i) {
			for (size_t j = 0; j + (size_t(1) << i) <= size(); ++j) {
				m_data[j][i] = m_oper(m_data[j][i - 1],
									  m_data[j + (size_t(1) << (i - 1))][i - 1]);
			}
		}
	}

	inline void clear() noexcept{ m_data.clear(); }

	inline size_t size() const noexcept { return m_data.size(); }
	inline bool empty() const noexcept { return m_data.empty(); }

	inline Elem query(size_t pos, size_t len) const {
		if (pos >= m_data.size()) {
			throw std::out_of_range("SparseTable::query: pos (which is " +
									std::to_string(pos) +
									") > size() (which is " +
									std::to_string(m_data.size()) + ')');
		}
		if (len == 0) {
			throw std::out_of_range("SparseTable::query: len == 0");
		}
		if (pos + len > size()) {
			len = size() - pos;
		}
		size_t log_len = floorLogn2(len);
		return m_oper(m_data[pos][log_len],
					  m_data[pos + len - (size_t(1) << log_len)][log_len]);
	}

protected:
	Oper m_oper;
	Table m_data; // m_data[i][j] maintains Oper(data[i, i + 2 ** j))

private:
};

[2025-07-31] FFT和NTT

FFT

思路来自Reducible的视频《The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?》，代码优化来自OI-Wiki。

快速傅里叶变换（Fast Fourier Transform）是多项式乘法的优化算法，可以在$O(n\log n)$时间内完成两个$n$次多项式之积的系数计算。

问题

为了下标的统一，我们设输入为两个$(n-1)$次多项式$f(x)=\sum_{i=0}^{n-1}a_ix^i,\ g(x)=\sum_{i=0}^{n-1}b_ix^i$。于是，他们的系数表示就是向量$\vec a=(a_i)_{i=0}^{n-1}$和$\vec b=(b_i)_{i=0}^{n-1}$。我们需要求出其积的系数，即求出$\vec c\in\R^n$满足$h(x)=\sum_{i=0}^{n-1}c_ix^i=(f\cdot g)(x)$。

朴素方法

直接展开两个多项式，得到其积$(f\cdot g)(x)=\sum_{i=0}^{2n-2}\left(\sum_{j=0}^{\min\{i,n-1\}}a_jb_{i-j}\right)x^i$。时间复杂度$O(n^2)$。

第一个优化：点表示

现在先只考虑一个多项式。

特别地，当$n=2$时，其图像是一条直线，可以由不重合的两点确定。一般地，任何一个$(n-1)$次多项式都可以这样$n$个点唯一确定，它们满足其中任意两点没有相同的横坐标。

证明

设$\vec a,\vec x,\vec y\in\R^n$满足$(\forall i\in\mathbb N\cap[0,n)\ \forall j\in\mathbb N\cap[0,n)\ (x_i=x_j\iff i=j)$。

补充定义$0^0=1$。记$X=\left(\left[x_i^j\right]_{j=0}^{n-1}\right)_{i=0}^{n-1}$，即

\[X= \begin{bmatrix} x_0^0 & x_0^1 & x_0^2 & \cdots & x_0^{n-1}\\ x_1^0 & x_1^1 & x_1^2 & \cdots & x_1^{n-1}\\ x_2^0 & x_2^1 & x_2^2 & \cdots & x_2^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ x_{n-1}^0 & x_{n-1}^1 & x_{n-1}^2 & \cdots & x_{n-1}^{n-1}\\ \end{bmatrix} \]
现在证明，存在唯一的$(n-1)$次多项式$f(x)=\sum_{i=0}^{n-1}a_ix^i$，使得$\forall i\in\mathbb N\cap[0,n)\ f(x_i)=y_i$。

我们将这$n$条方程写成矩阵：$\vec y=X\vec a$，即

\[\begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n-1} \end{bmatrix} = \begin{bmatrix} x_0^0 & x_0^1 & x_0^2 & \cdots & x_0^{n-1}\\ x_1^0 & x_1^1 & x_1^2 & \cdots & x_1^{n-1}\\ x_2^0 & x_2^1 & x_2^2 & \cdots & x_2^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ x_{n-1}^0 & x_{n-1}^1 & x_{n-1}^2 & \cdots & x_{n-1}^{n-1}\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n-1} \end{bmatrix} \]
由于$(\forall i\in\mathbb N\cap[0,n)\ \forall j\in\mathbb N\cap[0,n)\ (x_i=x_j\iff i=j)$，故$X$第二列互不相同，而第一列全为$1$，所以行之间不可能线性相关，$\mathrm{rk\ }X=n$，$X$可逆。

由此可得$\vec a=X^{-1}\vec y$。由逆矩阵的唯一性可知，此时$\vec a$唯一。

题外话：同时可以注意到，如果有两个$x$相等，则$X$必有两行完全相等，必不可逆。所以我们取的其实是充要条件。

用这样$n$个点表示多项式的方法称为多项式的点表示法。上述证明过程也是系数表示法和点表示法相互转换的过程。

第一个优化就是，将$f$和$g$转换为点表示，从而求出$h$的点表示，再利用逆矩阵得到$h$的系数表示。

具体地：

取$(2n-2)$个互不相同的横坐标$\vec x\in\R^{2n-2}$，分别代入得$f(\vec x)\triangleq(f(x_i))_{i=0}^{2n-2}\ ,g(\vec x)\triangleq(g(x_i))_{i=0}^{2n-2}$。
做Hadamard积，得$h(\vec x)\triangleq(h(x_i))_{i=0}^{2n-2}=f(\vec x)\odot g(\vec x)=(f(x_i)\cdot g(x_i))_{i=0}^{2n-2}$。
得到了$h(\vec x)$后，再利用$X=\left(\left[x_i^j\right]_{j=0}^{n-1}\right)_{i=0}^{n-1}$的逆矩阵，计算$h$的系数。

以上过程看似让计算变得复杂了，但却提供了大量的优化空间。2. 的时间复杂度为$O(n)$，已经足够小，优化的重点主要在1. 和3. 。

第二个优化：函数的奇偶性

$x$的奇数次幂是奇函数，偶数次幂是偶函数。无论是奇函数还是偶函数，我们都可以利用其对称性，求出$y$轴一侧的值而知道另一侧的值。所以，为了减少取点的次数，我们考虑总是取一对相反数，即$\pm x_0,\pm x_1,...,\pm x_{n/2-1}\ (x_i>0)$。为了能让$2\mid n$，将$f$补为偶数项（最高次为奇数次）。同时，我们对$f(x)$的项奇偶分组：

\[\begin{align*} f(x)&=\sum_{i=0}^{n-1}a_ix^i\\ &=\sum_{i=0}^{n/2-1}a_{2i}x^{2i}+\sum_{i=0}^{n/2-1}a_{2i+1}x^{2i+1}\\ &=\sum_{i=0}^{n/2-1}a_{2i}(x^i)^2+x\sum_{i=0}^{n/2-1}a_{2i+1}(x^i)^2 \end{align*} \]

设$f_\mathrm e(x)=\sum_{i=0}^{n/2-1}a_{2i}x^i,\ f_\mathrm o(x)=\sum_{i=0}^{n/2-1}a_{2i+1}x^i$。注意到了吗？这时，计算$f(x)$已经转化为了规模减半的两个子问题，计算$f_\mathrm e(x^2)$和$f_\mathrm o(x^2)$。

整个问题就变成：取$\pm x_0,\pm x_1,...,\pm x_{n/2-1}\ (x_i>0)$，计算

\[f(x_i)=f_\mathrm e(x_i^2)+x_if_\mathrm o(x_i^2)\\ f(-x_i)=f_\mathrm e(x_i^2)-x_if_\mathrm o(x_i^2) \]

总的时间复杂度减小为$O(n\log n)$。

不过，这里还有一个大问题。注意到了吗？$x_i^2>0$。这会导致子问题中我们无法继续分治。

x0  -x0  x1 -x1
  \  /    \  /
  x0**2  x1**2
     \   /
     x0**4

上图中，理想的情况是，x0**2与x1**2互为相反数，且$(x_0^2)^2=(x_1^2)^2=(x_0)^4$。显然，这在实数域内是不可能成立的，我们必须考虑

拓展到复数域

记$\omega_i\ (i\in\mathbb N^*)$为$i$次单位根，即$\omega_i=\exp\left(\frac{2\pi}{i}\mathrm i\right)$。这里我们要利用的性质是$\forall m\in\mathbb N^*\ (\forall i\in\mathbb N,\ i<2m)\ \omega_{2m}^i+\omega_{2m}^{i+m}=0$。

所以，我们要改一改上面的说法：不取$\pm x_0,\pm x_1,...,\pm x_{n/2-1}\ (x_i>0)$，而是取$\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$。于是有

\[f(\omega_n^i)=f_\mathrm e(\omega_n^{2i})+\omega_n^if_\mathrm o(\omega_n^{2i})=f_\mathrm e(\omega_{n/2}^i)+\omega_n^if_\mathrm o(\omega_{n/2}^i)\\ f(\omega_n^{i+n/2})=f_\mathrm e(\omega_n^{2i+n})+\omega_n^{i+n/2}f_\mathrm o(\omega_n^{2i+n})=f_\mathrm e(\omega_{n/2}^i)-\omega_n^if_\mathrm o(\omega_{n/2}^i) \]

现在还有一个小小的问题：必须让递归过程中，$n$要么是偶数，要么是$1$。所以一开始的时候，不仅要把$f$补为偶数项，而且要补到$2$的幂项，即要使$\log_2n\in\mathbb N$。

FFT

终于，我们真正将1. 优化到了$O(n\log n)$。

std::vector<Complex> fft(const std::vector<Complex> &coefs) { // assuming coefs.size() is a power of 2
	if (coefs.size() <= 1) {
		return coefs;
	}
	Complex omega = exp(Complex(0, 2 * M_PI / coefs.size()));
	auto even = std::vector<Complex>(coefs.size() >> 1),
		 odd = std::vector<Complex>(coefs.size() >> 1);
	for (size_t i = 0; (i << 1) != coefs.size(); ++i) {
		even[i] = coefs[i << 1];
		odd[i] = coefs[(i << 1) | 1];
	}
	even = fft(even);
	odd = fft(odd);
	auto res = std::vector<Complex>(coefs.size());
	Complex pow_omega = 1;
	for (size_t i = 0; (i << 1) != coefs.size(); ++i) {
		res[i] = even[i] + pow_omega * odd[i];
		res[i + (coefs.size() >> 1)] = even[i] - pow_omega * odd[i];
		pow_omega *= omega;
	}
	return res;
}

IFFT

现在我们需要将3. 也优化到$O(n\log n)$。

FFT的过程是计算$\vec y=X\vec a$，而现在我们需要计算$\vec a=X^{-1}\vec y$。这个过程称为逆快速傅里叶变换（IFFT）。

$X$完美符合Vandermonde矩阵的形式。这种矩阵有一些美妙的性质，例如我们熟悉的$\det X=\prod_{i=0}^{n-1}\prod_{j=i+1}^{n-1}(x_j-x_i)$。不过这里我们不需要复杂地推导一般的Vandermonde矩阵之逆，我们仅考虑$x_i=\omega_n^i$这种特殊情况。下面为了避免符号冲突，将代入了值的$X$记作$\Omega$。

\[\Omega= \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \cdots & (\omega_n^0)^{n-1}\\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1)^{n-1}\\ (\omega_n^2)^0 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2)^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1})^{n-1}\\ \end{bmatrix} \]

这个矩阵通常被称作离散傅里叶变换（DFT）矩阵。

可以求得

\[\Omega^{-1}=\frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & (\omega_n^{-0})^2 & \cdots & (\omega_n^{-0})^{n-1}\\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & (\omega_n^{-1})^2 & \cdots & (\omega_n^{-1})^{n-1}\\ (\omega_n^{-2})^0 & (\omega_n^{-2})^1 & (\omega_n^{-2})^2 & \cdots & (\omega_n^{-2})^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & (\omega_n^{-(n-1)})^2 & \cdots & (\omega_n^{-(n-1)})^{n-1}\\ \end{bmatrix} \]

推导

即使直接跳过这一部分也完全不会影响严谨性，这里纯粹只是满足一些好奇心。

共轭转置恰为逆矩阵的复方阵称为酉矩阵（unitary matrix）。严谨地，若矩阵$M\in\mathbb C^{n\times n}\ (n\in\mathbb N^*)$满足$\overline M^\top=M^{-1}$，则称$M$为酉矩阵。因$\Omega$是复矩阵，故考虑$\Omega$是否为酉矩阵。

\[\overline\Omega^\top= \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & (\omega_n^{-0})^2 & \cdots & (\omega_n^{-0})^{n-1}\\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & (\omega_n^{-1})^2 & \cdots & (\omega_n^{-1})^{n-1}\\ (\omega_n^{-2})^0 & (\omega_n^{-2})^1 & (\omega_n^{-2})^2 & \cdots & (\omega_n^{-2})^{n-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & (\omega_n^{-(n-1)})^2 & \cdots & (\omega_n^{-(n-1)})^{n-1}\\ \end{bmatrix}\\ \overline\Omega^\top\Omega=nI_n \]
所以就得到$\Omega^{-1}=\frac{1}{n}\overline\Omega^\top$。

另一种思路是，回到信号处理中，根据三角函数的正交性，构造推导逆离散傅里叶变换公式。由于涉及信号处理原理，这里就不不懂装懂了。

如果读者确实没有看上面一部分的推导，由于$\Omega$是满秩矩阵，我们也可以直接验证

\[\Omega^{-1}\Omega=I_n \]

验证

对任意$i,j\in\mathbb N\cap[0,n)$有

\[\begin{align*} (\Omega^{-1}\Omega)(i,j)&=\sum_{k=0}^{n-1}\Omega^{-1}(i,k)\Omega(k,j)\\ &=\sum_{k=0}^{n-1}\left(\frac{1}{n}\omega_n^{-i\cdot k}\right)\omega_n^{k\cdot j}\\ &=\sum_{k=0}^{n-1}\frac{1}{n}\omega_n^{k(j-i)} \end{align*} \]
同理于三角函数的正交性，复指数函数也有正交性：

\[\sum_{k=0}^{n-1}\exp\left(\frac{2\pi\mathrm i}{n}(j-i)k\right)= \left\{ \begin{align*} &0&(i\neq j)\\ &n&(i=j) \end{align*} \right. \ (i,j\in\mathbb Z) \]

证明

若$i=j$，则$原式=\sum_{k=0}^{n-1}1=n$。

若$i\neq j$，设$q=\exp\left(\frac{2\pi\mathrm i}{n}(j-i)\right)$，则由等比数列求和

\[\begin{align*} 原式&=\sum_{k=0}^{n-1}q^k\\ &=\frac{q^n-1}{q-1} \end{align*} \]
因$q^n=\exp(2\pi\mathrm i(j-i))$，其中$j-i\in\mathbb Z^*$，则由Euler公式可知$q^n=\cos(2\pi(j-i))+\sin(2\pi(j-i))=1$，故$原式=0$。

利用复指数函数的正交性即得$\Omega^{-1}\Omega=I_n$。

注意到$\Omega^{-1}$形式与$\Omega$非常相似，唯将$\omega_n$改为$\omega_n^{-1}$，并乘上归一化系数$\frac{1}{n}$。我们完全可以照搬FFT代码。

std::vector<Complex> ifft(const std::vector<Complex> &vals) { // assuming coefs.size() is a power of 2
	if (vals.size() <= 1) {
		return vals;
	}
	Complex omega = exp(Complex(0, -2 * M_PI / vals.size()));
	auto even = std::vector<Complex>(vals.size() >> 1),
		 odd = std::vector<Complex>(vals.size() >> 1);
	for (size_t i = 0; (i << 1) != vals.size(); ++i) {
		even[i] = vals[i << 1];
		odd[i] = vals[(i << 1) | 1];
	}
	even = ifft(even);
	odd = ifft(odd);
	auto res = std::vector<Complex>(vals.size());
	Complex pow_omega = 1;
	for (size_t i = 0; (i << 1) != vals.size(); ++i) {
		res[i] = even[i] + pow_omega * odd[i];
		res[i + (vals.size() >> 1)] = even[i] - pow_omega * odd[i];
		pow_omega *= omega;
	}
	return res;
}

注意返回值是没有归一化的。

递归实现FFT/IFFT

将“是否为逆运算”作为参数传入函数，可以大幅减少代码重复。

std::vector<Complex> fft(const std::vector<Complex> &arr, bool inv = false) { // assuming arr.size() is a power of 2
	if (arr.size() <= 1) {
		return arr;
	}
	Complex omega = exp(Complex(0, (inv ? -2 : 2) * M_PI / arr.size()));
	auto even = std::vector<Complex>(arr.size() >> 1),
		 odd = std::vector<Complex>(arr.size() >> 1);
	for (size_t i = 0; (i << 1) != arr.size(); ++i) {
		even[i] = arr[i << 1];
		odd[i] = arr[(i << 1) | 1];
	}
	even = fft(even, inv);
	odd = fft(odd, inv);
	auto res = std::vector<Complex>(arr.size());
	Complex pow_omega = 1;
	for (size_t i = 0; (i << 1) != arr.size(); ++i) {
		res[i] = even[i] + pow_omega * odd[i];
		res[i + (arr.size() >> 1)] = even[i] - pow_omega * odd[i];
		pow_omega *= omega;
	}
	return res;
}

注意执行IFFT时，返回值是没有归一化的。

位逆序置换优化

递归计算会占用大量空间。我们希望找到一种办法免于递归。

以$8$项多项式为例，模拟拆分的过程：

初始序列为$\{x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7\}$。
一次二分之后：$\{x_0,x_2,x_4,x_6\},\{x_1,x_3,x_5,x_7\}$。
两次二分之后：$\{x_0,x_4\}\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\}$。
三次二分之后$\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7\}$。

规律：其实就是原来的那个序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。比如$x_1$是$001$，翻转是$100$，也就是$4$，而且最后那个位置确实是$4$。我们称这个变换为位逆序置换（bit-reversal permutation)。

证明

第$i$次二分，就是在每个已经分好的块内按照从低到高第$i$个二进制位升序排序。这意味着，总的顺序定义为，从低位到高位比较，相同则跳过，否则较小的在前面。因此，最终的顺序就是下标的二进制表示翻转后从小到大的顺序。

朴素的方法是直接暴力枚举、翻转，时间复杂度$O(n\log n)$。

设$r(n,x)$为$n$位二进制数$x$翻转后的数。利用位运算，可以在$O(2^n)$时间内求出所有$r(n,x)\ (x\in\mathbb N\cap[0,2^n)$。递推公式是r(n,x) = (r(n, x >> 1) >> 1) | ((x & 1) << (n - 1))。

template <typename Iter>
using RequireFwdIter = std::enable_if_t<
	std::is_base_of<
		std::forward_iterator_tag,
		typename std::iterator_traits<Iter>::iterator_category>::value>;

template <typename Iter, typename = RequireFwdIter<Iter>>
inline void bitRevPerm(Iter begin, Iter end) {
	size_t n = std::distance(begin, end);
	if (n & (n - 1)) {
		throw std::invalid_argument("bitRevPerm: std::distance(begin, end)"
									" is not a power of 2");
	}
	auto rev = std::vector<size_t>(n);
	for (size_t i = 1; i < n; ++i) {
		rev[i] = (rev[i >> 1] >> 1);
		if (i & 1) {
			rev[i] |= (n >> 1);
		}
	}
	Iter iter = begin;
	for (size_t i = 0; i != n; ++i, ++iter) {
		if (i < rev[i]) {
			std::swap(*iter, *std::next(begin, rev[i]));
		}
	}
}

蝶形运算优化

现在我们想办法重复利用空间，而不是每次都重新开两个数组。

注意我们上面实现的fft函数push up的部分：

res[i] = even[i] + pow_omega * odd[i];
res[i + (arr.size() >> 1)] = even[i] - pow_omega * odd[i];

注意到了么？第$i$项和第$(i+n/2)$项的计算只涉及两个子问题结果的第$i$项，即$f(\omega_n^i)$和$f(\omega_n^{i+n/2})$都只与$f_\mathrm e(\omega_{n/2}^i)$和$f_\mathrm o(\omega_{n/2}^i)$有关。这意味着我们可能可以“原地修改”。

设当前$f$保存到的区间为$[0,n)$。我们希望$f_\mathrm e$和$f_\mathrm o$的结果分别占据$[0,n/2)$和$[n/2,n)$，因为这样很方便递归上下两层之间交互；另一方面，$f$的结果也应当按顺序保存在$[0,n)$，以便上一层访问。当$\left(f_\mathrm e(\omega_{n/2}^i)\right)_{i=0}^{n/2-1}$依次存储在$[0,n/2)$，$\left(f_\mathrm o(\omega_{n/2}^i)\right)_{i=0}^{n/2-1}$依次存储在$[n/2,n)$时，我们需要提取位置$i$和$n/2+i$上的值，计算后再分别放入这两个位置。多么和谐！

非递归实现FFT

inline void fft(std::vector<Complex> &arr, bool inv = false) {
	bitRevPerm(arr.begin(), arr.end());
	for (size_t len = 2; len <= arr.size(); len <<= 1) {
		Complex omega = exp(Complex(0, (inv ? -2 : 2) * M_PI / len));
		for (size_t i = 0; i != arr.size(); i += len) {
			Complex omega_pow = 1;
			for (size_t j = 0; (j << 1) != len; ++j) {
				Complex even = arr[i + j],
						odd = omega_pow * arr[i + (len >> 1) + j];
				arr[i + j] = even + odd;
				arr[i + (len >> 1) + j] = even - odd;
				omega_pow *= omega;
			}
		}
	}
	if (inv) {
		for (auto &x : arr) {
			x /= arr.size();
		}
	}
}

注意，这里的返回值无需再做归一化处理。

时间复杂度：$O(n\log n)$。
空间复杂度：原地修改，$O(1)$额外空间。

终于结束了！总共三次试图学习FFT，估计总用时近10小时。我真是太菜了。

NTT

费马小定理

Fermat's Little Theorem: 设$p$为质数，则$\forall a\in\mathbb Z\ (a^p\equiv a\pmod p\wedge(p\nmid a\implies a^{p-1}\equiv 1\pmod p))$。

证明

可以用数学归纳法证明。

当$a=0$时，$a^p=0\equiv a\pmod p$。

假设$a^p\equiv a\pmod p$，则

\[(1+a)^p=\sum_{i=0}^{p}\binom{p}{i}a^i=\sum_{i=0}^p\frac{\prod_{j=0}^{i-1}(p-j)}{i!}a^i\equiv1+a^p\equiv1+a\pmod p \]
故$\forall a\in\mathbb N\ (a^p\equiv a\pmod p)$。

再来考虑$\mathbb Z^-$的部分。

若$p=2$，由于平方不改变奇偶性，显然成立。下面考虑$p\neq 2$。

由于$p$是质数，所以$2\nmid p$，故$\forall i\in\mathbb Z\ ((p-i)\equiv i+1\pmod p)$。

假设$a^p\equiv a\pmod p$，则

\[(a-1)^p=\sum_{i=0}^p(-1)^{p-i}\binom{p}{i}a^i=\sum_{i=0}^p(-1)^{i+1}\frac{\prod_{j=0}^{i-1}(p-j)}{i!}a^i\equiv a^p-1\equiv a-1\pmod p \]
综上，$\forall a\in\mathbb Z\ (a^p\equiv a\pmod p)$。特别地，如果$p\nmid a$，则由$a(a^{p-1}-1)\equiv0\pmod p$可得$a^{p-1}\equiv1\pmod p$。

[2025-08-02] 数论基础

今天睡觉落枕了，脖子不可屈伸，又睡过头，怠了半天。

为了公式和代码之简便，我们做一些违背祖宗的定义：

定义$a/b$为与$\frac{a}{b}$相邻的整数中靠近$0$的那个。
定义$a\bmod b$为$a-a/b\cdot b$。
定义$\gcd(a,b)$为$a$和$b$的所有公约数中，与$b$符号相同且绝对值最大的那个。

乘法逆元

设$m\in\mathbb N^*,\ a\in\mathbb Z$，若$x\in\mathbb Z$满足$ax\equiv1\pmod m$，则称$x$为模$m$意义下$a$的乘法逆元。

扩展欧几里得算法

欧几里得算法（Euclidean algorithm）/辗转相除法：$\forall a,b\in\mathbb Z\ (\gcd(a,b)=\gcd(b,a\bmod b))$。

证明

设$a,b\in\mathbb Z,\ q=a/b,\ r=a-bq$。

任取$d\in\mathbb Z$满足$d\mid a\wedge d\mid b$（这样的$d$一定存在，例如$1$），则有$d\mid r$即$d\mid(a\bmod b)$。故$\{d\in\mathbb Z:\ d\mid a\wedge d\mid b\}\subseteq\{d\in\mathbb Z:\ d\mid b\wedge d\mid(a\bmod b)\}$。

任取$d\in\mathbb Z$满足$d\mid b\wedge d\mid(a\bmod b)$，即$d\mid b\wedge d\mid r$，故也有$d\mid(bq+r)$即$d\mid a$。故$\{d\in\mathbb Z:\ d\mid b\wedge d\mid(a\bmod b)\}\subseteq\{d\in\mathbb Z:\ d\mid a\wedge d\mid b\}$。

所以$\{d\in\mathbb Z:\ d\mid a\wedge d\mid b\}=\{d\in\mathbb Z:\ d\mid b\wedge d\mid(a\bmod b)\}$，$\gcd(a,b)=\max\{d\in\mathbb Z:\ d\mid a\wedge d\mid b\}=\max\{d\in\mathbb Z:\ d\mid b\wedge d\mid(a\bmod b)\}=\gcd(b,b\bmod a)$。正确性得证。

当$|a|<|b|$时，$(a,b)\gets(b,a\bmod b)$等价于$(a,b)\gets(b,a)$。

当$|a|=|b|$时，$(a,b)\gets(b,a\bmod b)$等价于$(a,b)\gets(b,0)$。此时$a$即为$\gcd(a,b)$。

当$|a|>|b|$时，$(a,b)\gets(b,a\bmod b)$会使$|a|,|b|$均减小。

综上，在有限次迭代后一定能得到$(a,b)=(a,0)$，此时$a$即为所求。可行性得证。

int64_t gcd(int64_t a, int64_t b) {
	if (!b) {
		return a;
	}
	return gcd(b, a % b);
} // 递归实现

int64_t gcd(int64_t a, int64_t b) {
	while (!b) {
		std::tie(a, b) = std::tuple(b, a % b);
	}
	return a;
} // 循环实现

定义$\mathrm{exGcd}(a,b)=(x,y)$，其中$x,y$满足$ax+by=\gcd(a,b)$。扩展欧几里得算法可以求出该方程的一组解。

扩展欧几里得算法（Extended Euclidean algorithm）：

$\mathrm{exGcd}(a,0)=(1,0)$。
设$(x,y)=\mathrm{exGcd}(b,a\bmod b)$，则$\mathrm{exGcd}(a,b)=(y,x-a/b\cdot y)$。

证明

当$b=0$时，显然可取$x=1,y=0$。

设$(x',y')=\mathrm{exGcd}(a,b),\ (x,y)=\mathrm{exGcd}(b,a\bmod b)$，即

\[ax'+by'=\gcd(a,b)\\ bx+(a\bmod b)y=\gcd(b,a\bmod b)=\gcd(a,b) \]
因此

\[ax'+by'=bx+(a\bmod b)y=bx+(a-a/b\cdot b)y \]
这里的除法向$0$舍入。移项，对$a,b$合并系数得

\[(x'-y)a+(y'-x+a/b\cdot y)b=0 \]
因此可取$x'=y,\ y'=x-a/b\cdot y$。

int64_t exGcd(int64_t a, int64_t b, int64_t &x, int64_t &y) {
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	auto d = exGcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
} // 递归实现

下面考虑将其改为循环实现。一次辗转相除，写成矩阵形式即为

\[\begin{bmatrix} b\\ a\bmod b \end{bmatrix} = \begin{bmatrix} & 1\\ 1 & -a/b \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} \]

如果设$M(a,b)=\begin{bmatrix} & 1\\1 & -a/b\end{bmatrix}$，则每次迭代就是

\[(a,b)\gets M(a,b)(a,b) \]

最终

\[(\gcd(a,b),0)=(M\circ M\circ\cdots\circ M)(a,b)(a,b)\triangleq \begin{bmatrix} x & y\\ u & v \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} \]

其中$x,y$即为所求。

那末，我们只需要维护这个系数矩阵的变化，即可在辗转相除的同时求出$x,y$。

初始时令$\begin{bmatrix}x & y\\u & v\end{bmatrix}=\begin{bmatrix}1 & \\ & 1\end{bmatrix}$即单位矩阵。每次迭代时，

\[\begin{bmatrix} x & y\\ u & v \end{bmatrix} \gets \begin{bmatrix} & 1\\ 1 & -a/b \end{bmatrix} \begin{bmatrix} x & y\\ u & v \end{bmatrix} \]

然后

\[(a,b)\gets(b,a\bmod b) \]

由于$a\bmod b=a-a/b\cdot b$，这里可以计算出$a/b$并重复利用。

int64_t exGcd(int64_t a, int64_t b, int64_t &x, int64_t &y) {
	x = 1, y = 0;
	int64_t u = 0, v = 1;
	while (b) {
		int64_t q = a / b;
		std::tie(a, b, x, y, u, v) =
			std::make_tuple(b, a - q * b, u, v, x - q * u, y - q * v);
	}
	return a;
} // 循环实现

时间复杂度：$O(\log\max\{a,b\})$。
空间复杂度：$O(1)$。

模$m$意义下$a$的乘法逆元存在的充要条件是$\gcd(a,m)=1$。

证明

充分性：求解$ax\equiv1\pmod m$只需解$ax+my=1=\gcd(a,m)$。利用exGcd可以直接解出$x$。

必要性：设$\gcd(a,m)\neq1$，则存在$d\in\mathbb Z\backslash\{1\}$满足$d\mid a\wedge d\mid m$。假设存在$x\in\mathbb Z$满足$ax\equiv1\pmod m$，则$ax=a/m\cdot m+1$。此式在模$d$意义下即为$0\equiv1\pmod d$，矛盾。所以假设不成立，此时$a$的乘法逆元不存在。

/**
 * @note assuming gcd(x, mod) == 1.
 */
inline int64_t modMulInv(int64_t x, int64_t mod) {
	int64_t res, tmp;
	exGcd(x, mod, res, tmp);
	return ((res % mod + mod) % mod);
}

时间复杂度：$O(\log\max\{a,m\})$。通常有$0\leq a<m$，此时为$O(\log m)$。
空间复杂度：$O(1)$。

快速幂

特别地，当$m$为质数时（此时一个整数只要不是$m$的倍数，它就一定有乘法逆元），由费马小定理，$\forall a\in\mathbb Z\ (m\nmid a\implies a^{m-1}\equiv1\pmod m)$。所以此时$a^{m-2}$即为$a$在模$m$意义下的乘法逆元。

时间复杂度：$O(\log m)$。
空间复杂度：$O(1)$。

递推法

仍然假设$m$为质数。如果需要求出所有$i\in\mathbb N\cap[1,n]$的逆元，直接应用上述算法的时间复杂度会达到$O(n\log m)$。

我们可以利用$m\bmod i<i$的性质实现递推。

当$i=1$时，$i^{-1}\equiv1\pmod m$。

现在我们假设$[1,i-1]$的逆元已经求出。由于$m\bmod i<i$，做带余除法分解$m=iq+r$。

\[m=iq+r\equiv0\pmod m \]

两边同时乘$i^{-1}r^{-1}$，得

\[qr^{-1}+i^{-1}\equiv0\pmod m \]

故

\[i^{-1}\equiv-qr^{-1}\pmod m \]

所以只需取

\[i^{-1}=(m-m/i)(m\bmod i)^{-1} \]

/**
 * @return modular multiplicative inverse of [0, @c max].
 * @note assuming @c mod is a prime; inverse of 0 is undefined.
 */
std::vector<uint32_t> modMulInvs(uint32_t max, uint64_t mod) {
	std::vector<uint32_t> invs(max + 1);
	invs[1] = 1;
	for (uint32_t i = 2; i <= max; ++i) {
		invs[i] = (mod - mod / i) * invs[mod % i] % mod;
	}
	return invs;
}

时间复杂度：$O(n)$。
空间复杂度：$O(1)$额外空间。

不过，上述做法只能求从$1$开始连续$n$个数的乘法逆元。要求任意$n$个数$\vec a$的乘法逆元，可以通过求前缀积$\vec p=\left(\left(\prod_{j=0}^{i-1}a_j\right)\bmod m\right)_{i=0}^{n}$。如果知道$p_{i+1}^{-1}\equiv\prod_{j=0}^{i}a_i^{-1}\pmod m$，则可以求得$a_i^{-1}\equiv a_i^{-1}\prod_{j=0}^{i-1}(a_j\cdot a_j^{-1})\equiv p_ip_{i+1}^{-1}\pmod m$，同时又可以求出$p_i^{-1}\equiv a_ip_{i+1}^{-1}\pmod m$。这个做法也不要求$m$为质数，只要$\forall i\in\mathbb N\cap[0,n)\ (\gcd(a_i,m)=1)$即可。

/**
 * @return modular multiplicative inverse of each element of @c arr.
 * @note assuming @c mod is a prime; inverse of 0 is undefined.
 */
std::vector<uint64_t> modMulInvs(const std::vector<uint64_t> &arr, uint64_t mod) {
	if (arr.empty()) {
		return arr;
	}
	auto pref_prods = std::vector<uint64_t>(arr.size() + 1);
	pref_prods[0] = 1;
	for (size_t i = 0; i != arr.size(); ++i) {
		pref_prods[i + 1] = pref_prods[i] * arr[i] % mod;
	}
	uint64_t prod_inv = modMulInv(pref_prods.back(), mod);
	auto invs = std::vector<uint64_t>(arr.size());
	for (size_t i = arr.size() - 1; ~i; --i) {
		invs[i] = pref_prods[i] * prod_inv % mod;
		prod_inv = prod_inv * arr[i] % mod;
	}
	return invs;
}

时间复杂度：$O(n+\log m)$。
空间复杂度：$O(n)$。

质数筛

埃氏筛（Sieve of Eratosthenes）

思路：对每个质数，筛去其倍数。

std::vector<uint32_t> sieveOfEratosthenes(uint32_t max) {
	if (max <= 1) {
		return {};
	}
	auto not_prime = std::vector<bool>(max + 1);
	std::vector<uint32_t> primes;
	primes.emplace_back(2);
	for (uint32_t i = 4; i <= max; i += 2) {
		not_prime[i] = true;
	}
	for (uint32_t i = 3; i * i <= max; i += 2) {
		if (not_prime[i]) {
			continue;
		}
		for (uint32_t j = (i << 1); j <= max; j += i) {
			not_prime[j] = true;
		}
	}
	for (uint32_t i = 3; i <= max; i += 2) {
		if (!not_prime[i]) {
			primes.emplace_back(i);
		}
	}
	return primes;
}

时间复杂度：$O(n\log\log n)$。
空间复杂度：$O(n)$。

欧氏/线性筛（Euler’s/Linear Sieve）

思路：对每个数，筛取其质数倍的数。

std::vector<uint32_t> linearSieve(uint32_t max) {
	if (max <= 1) {
		return {};
	}
	auto not_prime = std::vector<bool>(max + 1);
	std::vector<uint32_t> primes;
	primes.emplace_back(2);
	for (uint32_t i = 3; i <= max; i += 2) {
		if (!not_prime[i]) {
			primes.emplace_back(i);
		}
		for (auto j : primes) {
			if (j * i > max) {
				break;
			}
			not_prime[j * i] = true; // j is the minimum prime factor of i * j
			if (i % j == 0) {
				break;
			}
		}
	}
	return primes;
}

时间复杂度：$O(n)$。
空间复杂度：$O(n)$。

[2025-08-03] 狄利克雷卷积、莫反和杜教筛

数论基础

同余类和剩余系

为表达简便起见，对于集合$S,T\subseteq\mathbb C$和$x\in\mathbb C$，定义$x+S\triangleq\{x+s:s\in S\},\ xS\triangleq\{xs:s\in S\},\,\ S+T\triangleq\{s+t:s\in S,t\in T\},\ ST\triangleq\{st:\ s\in S,t\in T\}$。

设$m\in\mathbb Z^*$，将模$m$同余关系的等价类（即余数相同的数都等价）称为模$m$的同/剩余类，把包含$x$的同余类记作$x\bmod m$，把这个划分记作$\mathbb Z_m$。特别地，对于$x\in\mathbb Z$满足$\gcd(x,m)=1$的情况，称$x\bmod m$为模$m$的既约同/剩余类，把既约同余类的集合记作$\mathbb Z_m^*$。

对于$\vec a\in\mathbb Z^m$，若$(\forall x\in\mathbb Z)\ (\exist i\in\mathbb N\cap[0,m))\ (x\equiv a_i\pmod m\wedge(\forall j\in\mathbb N\cap[0,m)\backslash\{i\}\ (x\not\equiv a_j\pmod m))$，则称$\vec a$这$m$个整数为模$m$的(完全)剩余系。

完全剩余系有这样的性质：设$m\in\mathbb Z^*$，$Z_m$是模$m$的一个完全剩余系，则$\forall a\in\mathbb Z\ (\gcd(a,m)=1\implies aZ_m是模m的一个完全剩余系)$。

证明

可以通过抽屉原理证明。

不妨令$m>0$。设$\vec r\in\mathbb Z^m$满足$\forall i\in\mathbb N\cap[0,m)\ (r_i\in Z_m\ \wedge r_i\equiv i\pmod m)$。

先证明$|aZ_m|=m$。如果$a=0$，必有$\gcd(a,m)=m=1$，显然成立。当$a\neq0$时，由于$r_i$各不相同，故$ar_i$也各不相同，所以$|aZ_m|=m$。

再来证明$aZ_m$的元素在模$m$下互不同余。假设存在$i,j\in\mathbb N\cap[0,m)$满足$i\neq j\wedge a\cdot r_i\equiv a\cdot r_j\pmod m$，则$a(r_i-r_j)\equiv0\pmod m$。由于$\gcd(a,m)=1$，所以$r_i-r_j\equiv0\pmod m$即$r_i\equiv r_j\pmod m$，这与$Z_m$是模$m$的剩余系矛盾，故假设不成立，所以$aZ_m$的元素在模$m$下互不同余。

综上可知$aZ_m$是模$m$的完全剩余系。

对于同余类$r\bmod m$，若$\gcd(r,m)=1$，则该同余类的所有元素均与$m$互质，称该同余类为模$m$的既约同/剩余类。把模$m$的既约同余类的个数记为$\varphi(m)$，称为欧拉函数（Euler's totient function）。把所有模$m$的既约同余类的集合记作$\mathbb Z_m^*$。

对于$\vec a\in\mathbb Z^{\varphi(m)}$，若$\forall i\in\mathbb N\cap[0,n)\ (\gcd(a_i,m)=1)$，且$(\forall x\in\mathbb Z,\ x\bmod m=1)\ (\exist i\in\mathbb N\cap[0,m))\ (x\equiv a_i\pmod m\wedge(\forall j\in\mathbb N\cap[0,m)\backslash\{i\}\ (x\not\equiv a_j\pmod m))$，则称$\vec a$这$\varphi(m)$个整数为模$m$的既约/缩/简化剩余系。

剩余系的复合：设$m=\prod_{i=0}^{n-1}m_i$满足$\forall i\in\mathbb N\cap[0,n)\ (m_i\geq1)$，$Z_{m_i}$为模$m_i$的(完全)剩余系，则$Z_m=\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)Z_{m_i}$为模$m$的(完全)剩余系。

证明

仍然可以通过抽屉原理证明。从每个剩余系中选取一个元素，一共有$\prod_{i=0}^{n-1}m_i$种选取方式，这意味着$|Z_m|\leq\prod_{i=0}^{n-1}m_i$。如果任意两种选取方式产生的数都不会在模$m$意义下同余，该式即可取等，同时也恰好符合完全剩余系的定义。

设在不相同的两种选取方式下，选取的数分别为$\vec a,\vec b$，产生的数分别为$x,y$。更严谨地，设$\vec a,\vec b\in\mathbb Z^n$满足$\forall i\in\mathbb N\cap[0,n)\ (a_i\in Z_{m_i})$且$\vec a\neq \vec b$，设

\[x=\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)a_i\\ y=\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)b_i \]
下面用反证法证明$x\not\equiv y\pmod m$。假设$x\equiv y\pmod m$，即

\[\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)a_i\equiv \sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)b_i\pmod m \]
移项，合并系数得

\[\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)(a_i-b_i)\equiv0\pmod m\tag0 \]
这意味着

\[m\mid\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)(a_i-b_i) \]
如果一个数能被$m$整除，那这个数一定能被$m$的任何一个约数整除。由于$m=\prod_{i=0}^{n-1}m_i$，可以知道

\[\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)(a_i-b_i)\equiv0\pmod{m_0} \]
当$i\geq1$时，系数中都含有$m_0$，所以有

\[\sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)(a_i-b_i)\equiv a_0-b_0\equiv0\pmod{m_0} \]
因$a_0,b_0\in Z_{m_0}$，在完全剩余系$Z_{m_0}$中不会有不相等的两个数在模$m_0$意义下同余，故$a_0=b_0$。

于是$(0)$式化为

\[\sum_{i=1}^{n-1}\left(\prod_{j=0}^{i-1}m_j\right)(a_i-b_i)\equiv0\pmod{\left(\prod_{i=0}^{n-1}m_i\right)} \]
由于$m_0$同时是等式左边、等式右边和模数的因子，所以可以同时除以$m_0$得

\[\sum_{i=1}^{n-1}\left(\prod_{j=1}^{i-1}m_j\right)(a_i-b_i)\equiv0\pmod{\left(\prod_{i=1}^{n-1}m_i\right)}\tag1 \]
同理于$(0)$式可以得到$a_1=b_1$。重复上述过程就可以得到$\forall i\in\mathbb N\cap[0,n)\ (a_i=b_i)$即$\vec a=\vec b$，与假设矛盾，故假设不成立。

所以，两种选取得到的数在模$m$意义下一定不同余，故$|Z_m|=\prod_{i=0}^{n-1}m_i$，且$Z_m$中任意两元素在模$m$意义下不同余，所以$Z_m$是模$m$的完全剩余系。

既约剩余系的复合：设$m=m_0m_1$满足$m_0,m_1\geq1$且$\gcd(m_0,m_1)=1$，设$Z_{m_0},Z_{m_1}$分别为模$m_0,m_1$的既约剩余系，则$Z_m\triangleq m_1Z_{m_0}+m_0Z_{m_1}$为模$m$的既约剩余系。

证明

设$Z_{m_0}',Z_{m_1}'$分别为模$m_0,m_1$的完全剩余系。由于$\gcd(m_0,m_1)=1$，故$m_1Z_0$是模$m_0$的完全剩余系，故$Z_{m}'\triangleq m_1Z_{m_0}'+m_0Z_{m_1}'$是模$m$的完全剩余系。下面证明$\{x\in Z_m':\gcd(x,m)=1\}=Z_m$。将左边这个集合记作$X$。

因$Z_{m_0}\subseteq Z_{m_0}',\ Z_{m_1}\subseteq Z_{m_1}'$，可得$Z_m\subseteq Z_m'$。

取$(r_0,r_1)\in Z_{m_0}'\times Z_{m_1}'$满足$\gcd(m_1r_0+m_0r_1,m_0m_1)=1$，则$r\triangleq m_1r_0+m_0r_1\in X$。如果$r$与一个数互质，则$r$必与这个数的所有约数都互质。由$\gcd(r_0,r_1)=1$，应用欧几里得定理得

\[1=\gcd(r,m_0)=\gcd(m_0,r\bmod m_0)=\gcd(m_0,(m_1r_0)\bmod m_0)=\gcd(m_1r_0,m_0)=\gcd(r_0,m_0)\\ 1=\gcd(r,m_1)=\gcd(m_1,r\bmod m_1)=\gcd(m_1,(m_0r_1)\bmod m_1)=\gcd(m_0r_1,m_1)=\gcd(r_1,m_1) \]
于是知$(r_0,r_1)\in Z_{m_0}\times Z_{m_1}$，所以$r\in Z_m$。因为$r$是任取的，所以

\[X\subseteq Z_m \]
任取$(r_0,r_1)\in Z_{m_0}\times Z_{m_1}$，则$r\triangleq m_1r_0+m_0r_1\in Z_{m}$。由于

\[\gcd(r_0,m_0)=\gcd(r_1,m_1)=\gcd(m_0,m_1)=1 \]
所以

\[\gcd(m_1r_0,m_0)=\gcd(m_0r_1,m_1)=1 \]
进而

\[\gcd(r,m)=\gcd(m_1r_0+m_0r_1,m_0m_1)=1 \]
又$r\in Z_m,\ Z_m\subseteq Z_m'$，$r$是任取的，所以

\[Z_m\subseteq X \]
综上，$X=Z_m$。

数论函数

若函数$f:\ \mathbb N^*\to\mathbb C$，则称$f$为一个数论/算数函数。

积性函数

若数论函数$f$满足$f(1)=1,\ \forall x,y\in\mathbb N^*\ (\gcd(x,y)=1\implies f(xy)=f(x)f(y))$，则称$f$为积性函数（multiplicative function）。

特别地，若积性函数$f$满足$\forall x,y\in\mathbb N^*\ (f(xy)=f(x)f(y))$，则称$f$为完全积性函数。

若$f,g$均为积性函数，$p\in\R$，则以下函数也为积性函数：

\[h:x\mapsto f(x^p)\\ h=f^p\\ h=f\cdot g\\ h=f*g \]

定义几个常用的函数：

单位函数（完全积性）：$\varepsilon:n\mapsto(n=1)$。
幂函数（完全积性）：$I_p:n\mapsto n^p$。定义$I=I_1$。
常数函数（完全积性）：$1$。
约数函数（积性）：$\sigma_p:\sum_{d\in\mathbb N^*,d\mid n}d^p$。
欧拉函数（积性）：$\varphi:n\mapsto\sum_{i=1}^n(\gcd(i,n)=1)$。
莫比乌斯函数（积性）：$\mu:n\to\begin{cases}1&&(n=1)\\0&&(n含有平方因数)\\(-1)^{n的不相等质因数个数}&&(n做质因数分解后，各指数均为1)\end{cases}$。

狄利克雷卷积

对于数论函数$f,g$，称函数$h:x\mapsto\sum_{y,z\in\mathbb N^*,yz=x}f(y)\cdot g(z)$为$f$与$g$的迪利克雷卷积（Dirichlet convolution），记作$f*g$。

设$f,g,h$均为数论函数，则有运算律：

交换律：$f*g=g*f$。
结合律：$(f*g)*h=f*(g*h)$。

证明

设$n\in\mathbb N^*$。

\[\begin{align*} ((f*g)*h)(n)&=\sum_{x,y\in\mathbb N^*,\ xy=n}(f*g)(x)\cdot h(y)\\ &=\sum_{x,y\in\mathbb N^*,\ xy=n}h(y)\sum_{s,t\in\mathbb N^*,\ st=x}f(s)\cdot g(t)\\ &=\sum_{x,y,z\in\mathbb N^*,\ xyz=n}f(x)\cdot g(y)\cdot h(z)\\ &=\sum_{x,y\in\mathbb N^*,\ xy=n}f(x)\sum_{s,t\in\mathbb N^*,\ st=y}g(s)\cdot h(t)\\ &=\sum_{x,y\in\mathbb N^*,\ xy=n}f(x)\cdot(g*h)(y)\\ &=(f*(g*h))(n) \end{align*} \]
对加法的分配律：$f*(g+h)=f*g+f*h$。

证明

设$n\in\mathbb N^*$。

\[\begin{align*} (f*(g+h))(n)&=\sum_{x,y\in\mathbb N^*,\ xy=n}f(x)\cdot (g+h)(y)\\ &=\sum_{x,y\in\mathbb N^*,\ xy=n}f(x)\cdot g(y)+\sum_{x,y\in\mathbb N^*,\ xy=n}f(x)\cdot h(y)\\ &=(f*g+f*h)(n) \end{align*} \]
消去律：若$h(1)\neq0$，则$f*h=g*h\iff f=g$。

证明

充分性显然成立，下面证明必要性。

设$f*h=g*h$，则由分配律有

\[(f-g)*h=0 \]
假设$f\neq g$，令$n_0=\min\{n\in\mathbb N^*:f(n)\neq g(n)\}$，则有

\[((f-g)*h)(n_0)=\sum_{x,y\in\mathbb N^*,\ xy=n_0}(f-g)(x)\cdot h(y) \]
当$x< n_0$时，$(f-g)(x)=0$，所以

\[((f-g)*h)(n_0)=(f-g)(n_0)\cdot h(1)\neq 0 \]
与$(f-g)*h=0$矛盾，故假设不成立，$f=g$。
单位元：$\varepsilon*f=f$。
逆元：TODO

莫比乌斯反演

TODO

欧拉函数

欧拉函数是积性函数。

证明

在“同余类和剩余系”中提到既约剩余系的复合：设$m=m_0m_1$满足$m_0,m_1\geq1$且$\gcd(m_0,m_1)=1$，设$Z_{m_0},Z_{m_1}$分别为模$m_0,m_1$的既约剩余系，则$Z_m\triangleq m_1Z_{m_0}+m_0Z_{m_1}$为模$m$的既约剩余系。

在这个过程中，有$|Z_m|=|m_1Z_{m_0}+m_0Z_{m_1}|=|Z_{m_0}||Z_{m_1}|$。欧拉函数$\varphi(m)$就是模$m$的既约剩余系的大小（基数），由此可见$\varphi(m)=\varphi(m_0)\varphi(m_1)$。

[2025-08-04] BITACM SC 2025 Rating #3

组合数学

错排公式

设$d_n$表示$\mathbb N\cap[1,n]$全部错位的排列方案数，则有递推式

\[d_n=d_{n-1}n+(-1)^n\ (n\in\mathrm N^*) \]

递推边界

\[d_0=1 \]

有表达式

\[d_n=\left\lceil\frac{n!}{\mathrm{e}}\right\rceil-n\bmod 2\ (n\in\mathbb N) \]

inline uint64_t derangement(uint64_t n, uint64_t mod) {
	uint64_t d = 1;
	for (uint64_t i = 1; i <= n; ++i) {
		d = (i * d % mod + ((i & 1) ? (mod - 1) : 1)) % mod;
	}
	return d;
}

inline std::vector<uint64_t> derangements(uint64_t max, uint64_t mod) {
	auto d = std::vector<uint64_t>(max + 1);
	d[0] = 1;
	for (uint64_t i = 1; i <= max; ++i) {
		d[i] = (i * d[i - 1] % mod + ((i & 1) ? (mod - 1) : 1)) % mod;
	}
	return d;
}

时间复杂度：$O(n)$。
空间复杂度：$O(1)$额外空间。

排列数和组合数

模意义下，排列数和组合数可以通过预处理阶乘和阶乘逆来快速计算。阶乘逆一方面可以通过递推法求出$\mathbb N\cap[1,n]$的逆之后求前缀积，另一方面也可以求出$n!$后逆向递推：

\[n!^{-1}\equiv(n+1)(n+1)!^{-1}\pmod m \]

快速取模

神秘的卡常奇技淫巧，称为Barrett reduce。

这个算法可以避免耗时较大的除法运算，适合需要多次对同一个模数取模的场景。总的思路是，将$\frac{1}{m}$乘上一个很大的数后取整，从而以较高的精度保存这个数，将除以$m$转换为乘$\frac{1}{m}$。

在$b\in\mathbb N\cap[2,+\infty)$进制下这个算法的步骤如下：

预先给定模数$m\in\mathbb N^*$，设$l=\lfloor\log_bm\rfloor+1$（即$m$在$b$进制下的位数），预处理出

\[\mu=b^{2l}/m \]
输入$x\in\mathbb N\cap[0,m^2)$。
计算近似商

\[q_0=x/b^l \]
计算更精确的近似商

\[q=\mu\cdot q_0/b^l \]
计算余数

\[r=x-mq \]
重复$r\gets r-m$直至$0\leq r<m$，得到$r=x\bmod m$。

证明

为区分真实值与估计值，设

\[q=x/m\\ r=x\bmod m\\ \tilde {q_0}=x/b^l\\ \tilde q=\mu\cdot\tilde{q_0}/b^l\\ \tilde r=x-m\tilde q \]
由整除的性质知道

\[\mu\leq\frac{b^{2l}}{m}\\ \tilde{q_0}\leq\frac{x}{b^l} \]
代入$\tilde q$得

\[\tilde q\leq\frac{b^{2l}}{m}\cdot\frac{x}{b^l}/b^l=\frac{x}{m}\cdot b^l/b^l=x/m=q \]
另一方面，

\[\mu>\frac{b^{2l}}{m}-1\\ \tilde{q_0}>\frac{x}{b^l}-1 \]
两式相乘，

\[\mu\cdot\tilde{q_0}>b^l\frac{x}{m}-\frac{b^{2l}}{m}-\frac{x}{b^l}+1 \]
两边同时除以$b^l$（因为$\tilde q$的计算需要）得

\[\frac{\mu\cdot\tilde{q_0}}{b^l}>\frac{x}{m}-\frac{b^l}{m}-\frac{x}{b^{2l}}+\frac{1}{b^l} \]
所以

\[\tilde q>\frac{\mu\cdot\tilde{q_0}}{b^l}-1>\frac{x}{m}-\left(\frac{b^l}{m}+\frac{x}{b^{2l}}\right)+\frac{1}{b^l}-1 \]
使用放缩法：

由$l=\lfloor\log_bm\rfloor+1$得$b^{l-1}\leq m<b^l$，因此$\frac{b^l}{m}\leq b$。

由$x<m^2,\ m<b^l$得$x<b^{2l}$，因此$\frac{x}{b^{2l}}<1$。

综上得

\[\frac{b^l}{m}+\frac{x}{b^{2l}}<b+1 \]
代回得

\[\tilde q>\frac{x}{m}-(b+1)+\frac{1}{b^l}-1>q-(b+2) \]
因为$b,q,\tilde q\in\mathbb N$，所以

\[q-(b+1)\leq\tilde q\leq q \]
进而

\[r=x-mq\leq\tilde r\leq x-mq+m(b+1)=r+m(b+1) \]
所以最多只需要做$(b+1)$次减法即可将$\tilde r$减到$[0,m)$范围。不过这一过程的正确性仍需要证明。由于减模数$m$不改变模$m$下的值，所以只需证明$\tilde r\equiv r\pmod m$。

观察

\[r=x-mq\\ \tilde r=x-m\tilde q \]
因为$q,\tilde q\in\mathbb N$，所以显然有

\[r\equiv x\equiv \tilde r\pmod m \]

AtCoder库的实现则进一步优化，扩大了$x$允许的范围，并省掉了$q_0$的计算：

预先给定模数$m\in\mathbb N^*$，设$l=\lfloor\log_bm\rfloor+1$（即$m$在$b$进制下的位数），预处理出

\[\mu=\left\lceil\frac{b^{2l}}{m}\right\rceil \]
输入$x\in\mathbb N\cap[0,b^{2l})$。
计算近似商

\[q=\mu\cdot x/b^{2l} \]
计算余数

\[r=x-mq \]
若$x<mq$则令

\[r\gets r+m \]

证明

符号约定与上一个证明相同。仍然通过放缩法确定估计值的范围。

\[\frac{\mu\cdot x}{b^{2l}}<\frac{\left(\frac{b^{2l}}{m}+1\right)x}{b^{2l}}=\frac{x}{m}+\frac{x}{b^{2l}}<q+1 \]
两边同时向下取整时，小于号应变为小于等于号，因为左右的值可能夹在相邻两整数之间：

\[\tilde q=\left\lfloor\frac{\mu\cdot x}{b^{2l}}\right\rfloor\leq\lfloor q+1\rfloor=q+1 \]
另一方面，

\[\tilde q>\frac{\mu\cdot x}{b^{2l}}-1\geq\frac{\frac{b^{2l}}{m}\cdot x}{b^{2l}}-1=\frac{x}{m}-1=q-1 \]
因为$q,\tilde q\in\mathbb N$，所以

\[\tilde q\geq q \]
综上可得

\[q\leq\tilde q\leq q+1 \]
所以

\[r-m=x-mq-m\leq\tilde r=x-m\tilde q\leq x-mq=r \]

namespace std {
	template <>
	struct is_unsigned<__uint128_t> : public true_type {};
}

template <typename Mod = uint32_t,
		  typename Mu = uint64_t,
		  typename Aux = __uint128_t,
		  typename = std::enable_if_t<std::is_unsigned<Mod>::value &&
									  std::is_unsigned<Mu>::value &&
									  std::is_unsigned<Aux>::value &&
									  sizeof(Mod) * 2 == sizeof(Mu) &&
									  sizeof(Mu) * 2 == sizeof(Aux)>>
class Barrett {
public:
	explicit Barrett(Mod mod) : m_mod(mod), m_mu((Mu)(-1) / mod + 1) {}

	void assign(Mod mod) {
		m_mod = mod;
		m_mu = (Mu)(-1) / mod + 1;
	}

	Mod umod() const { return m_mod; }

	Mod operator()(Mu x) const {
		Mu p = ((static_cast<Aux>(x) * m_mu) >> (sizeof(Mu) * 8)) * m_mod;
		return ((x < p) ? (x - p + m_mod) : (x - p));
	}
	Mod operator()(Mod lhs, Mod rhs) const {
		return operator()(static_cast<Mu>(lhs) * rhs);
	}

protected:
	Mod m_mod;
	Mu m_mu;

private:
}

[2025-08-06] KMP、字符串哈希和Manacher

前缀函数

KMP算法是一种能在线性时间内快速完成字符串匹配的算法，其关键在于在线性时间内计算前缀函数。

前缀函数的计算

设字符集为$\Sigma$，字符串$\vec s\in\Sigma^n\ (n\in\mathbb N)$。

若$l\in\mathbb N\cap[0,n]$，则称$s_{[0,l)}$为$\vec s$的长为$l$的前缀，$s_{[n-l,n)}$为$\vec s$的长为$l$的后缀。若$l\in\mathbb N\cap[0,n)$，则称$s_{[0,l)}$为$\vec s$的长为$l$的真前缀，$s_{[n-l,n)}$为$\vec s$的长为$l$的真后缀。若$\vec s$的长为$l$的真前缀和真后缀相等，则称$s_{[0,l)}$为$\vec s$的一个border。

定义$\vec s$的前缀函数为其各前缀的border长度，即

\[\vec\varpi(\vec s)\triangleq\left(\max_{l\in\mathbb N,\ l<n,\ s_{[0,l)}=s_{(i-l,i]}}l\right)_{i=0}^{n-1} \]

计算前缀函数的朴素方法需要三重循环，时间复杂度达到$O(n^3)$。这个过程中会有很多重复的比较，所以考虑通过递推来避免重复。

第一个优化

从$(i-1)$到$i$时，可以先只比较$s_{\varpi(\vec s)_{i-1}}$与新加进来的字符$s_i$，如果不匹配再重新算。

\[\varpi(\vec s)_i= \begin{cases} \varpi(\vec s)_{i-1}+1&(s_{\varpi(\vec s)_{i-1}}=s_i)\\ \max_{l\in\mathbb N,\ s_{[0,l)}=s_{(i-l,i]}}l&(s_{\varpi(\vec s)_{i-1}}\neq s_i) \end{cases} \]

第二个优化

对于新加的字符不匹配的情况，可以考虑寻找一个小于$\varpi(\vec s)_{i-1}$的最长的长度$j$满足$s_{[0,j)}=s_{(i-1-j,i-1]}$，我们再比较$s_j$与$s_i$。

因为

\[s_{[0,\varpi(\vec s)_{i-1})}=s_{(i-1-\varpi(\vec s)_{i-1},i-1]}\\ s_{[0,j)}=s_{(i-1-j,i-1]}\\ j<\varpi(\vec s)_{i-1} \]

所以

\[s_{[0,j)}=s_{[\varpi(\vec s)_{i-1}-j,\varpi(\vec s)_{i-1})} \]

我们希望$j$在小于$\varpi(\vec s)_{i-1}$的条件下越大越好，也就是说

\[j=\max_{l\in\mathbb N,\ l<\varpi(\vec s)_{i-1},\ s_{[0,l)}=s_{(\varpi(\vec s)_{i-1}-1-l,\varpi(\vec s)_{i-1}-1]}}l=\varpi(\vec s)_{\varpi(\vec s)_{i-1}-1} \]

那如果这个$j$也匹配不上，我们就需要取小于$j$的最大的$j'$使得$s_{[0,j')}=s_{(i-1-j',i-1]}$，如此重复。

\[\begin{array}{ll} 1 & j\gets\varpi(\vec s)_{i-1}\\ 2 & \textbf{while }j>0\text{ and }s_j\neq s_i\\ 3 & \qquad j\gets\varpi(\vec s)_{j-1}\\ 4 & \varpi(\vec s)_i\gets j+(s_j=s_i) \end{array} \]

template <typename Iter, typename = RequireFwdIter<Iter>>
inline std::vector<size_t> prefFuncOf(Iter begin, Iter end) {
	auto pi = std::vector<size_t>(std::distance(begin, end));
	end = std::next(begin);
	for (size_t i = 1, j; i < pi.size(); ++i, ++end) {
		for (j = pi[i - 1]; j && *end != *std::next(begin, j); j = pi[j - 1])
			;
		pi[i] = j + (*end == *std::next(begin, j));
	}
	return pi;
}

时间复杂度：$O(n)$。
空间复杂度：$O(1)$额外空间。

KMP算法

设$\Sigma$为字符集，给定$\vec s\in\Sigma^m,\ \vec t\in\Sigma^n\ (m,n\in\mathbb N)$，需要找出$s$在$t$中出现的位置，即求所有的$i\in\mathbb N\cap[0,n-m]$满足$t_{[i,i+m)}=\vec s$。

只需将$\vec s,\vec t$拼接，然后计算前缀函数即可。具体来说，设$\vec r=(\vec s,\mathtt{'\#'},\vec t)$，其中分隔符$\mathtt{'\#'}\not\in\Sigma$，求出$\varpi(\vec r)$，然后找到所有$i\in\mathbb N\cap[0,n-m]$满足$\varpi(\vec r)_{i+2m}=m$。

template <typename Iter>
inline std::vector<size_t>
kmp(Iter text_begin, Iter text_end,
	Iter pattern_begin, Iter pattern_end,
	const typename std::iterator_traits<Iter>::value_type &sep = -1) {
	size_t text_sz = std::distance(text_begin, text_end),
		   pattern_sz = std::distance(pattern_begin, pattern_end);
	auto seq = std::vector<typename std::iterator_traits<Iter>::value_type>(
		text_sz + 1 + pattern_sz);
	std::copy(pattern_begin, pattern_end, seq.begin());
	std::copy(text_begin, text_end, seq.end() - text_sz);
	seq[pattern_sz] = sep;
	auto pi = prefFuncOf(seq.begin(), seq.end());
	std::vector<size_t> res;
	for (size_t i = (pattern_sz << 1); i != pi.size(); ++i) {
		if (pi[i] == pattern_sz) {
			res.emplace_back(i - (pattern_sz << 1));
		}
	}
	return res;
}

时间复杂度：$O(m+n)$。
空间复杂度：$O(m+n)$。

计算序列的最小正周期

设$\Sigma$为字符集，字符串$\vec s\in\Sigma^n\ (n\in\mathbb N)$。

若$p\in\mathbb N^*\cap[1,n]$满足$\forall i\in\mathbb N\cap[0,n-p)\ (s_i=s_{i+p})$（最后一个循环可能不完全），则称$p$为$\vec s$的一个周期。

若$\vec s$有长为$l$的border，则$(n-l)$必为$\vec s$的一个周期。要求最小正周期，即求最大border长度，即求前缀函数。

template <typename Iter>
inline std::vector<size_t> borderLengths(Iter begin, Iter end) {
	auto pi = prefFuncOf(begin, end);
	std::vector<size_t> res;
	for (size_t i = pi.size(); i; i = pi[i - 1]) {
		res.emplace_back(pi[i - 1]);
	}
	return res;
}

/**
 * @note the last period may be incomplete
 */
template <typename Iter>
inline size_t periodLength(Iter begin, Iter end) {
	if (begin == end) {
		return 0;
	}
	auto pi = prefFuncOf(begin, end);
	return (pi.size() - pi.back());
}

时间复杂度：$O(n)$。
空间复杂度：$O(n)$。

如果要求$p\mid n$，则只需((s.size() % p) ? s.size() : p)即可。

证明

引理：周期引理（Fine and Wilf's theorem/periodicity lemma）：设Σ为字符集，字符串$\vec s\in\Sigma^n\ (n\in\mathbb N)$。若$p,q$都是$\vec s$的周期，$n\geq(p+q)-\gcd(p,q)$，则$\gcd(p,q)$也是$\vec s$的周期。

这个定理在1963年才被提出？这实在令人惊奇，我以为古代数学家肯定早就证明过，因为他的形式实在非常优美简洁；还是说，其实提出过，只是没有被记载？

证明

没有找到一个我看得懂的证明。不过有一个弱周期引理（weak periodicity lemma）很容易证明，只需将上述定理中$n$的范围改成$n\geq p+q$即可。

GPT-5.0的推出为我带来了更加逻辑清晰且简洁的证明，尽管其仍有思维混乱之处，但为我指明了大方向。下面我们就来证明周期引理。总体的思路是按照更相减损术的思路递归证明。

设命题

\[cond_0(b,e,p,q):(e-b\geq(p+q)-\gcd(p,q))\\ cond_1(b,e,p,q):p,q均为s_{[b,e)}的周期\\ cond(b,e,p,q):(cond_0(b,e,p,q)\wedge cond_1(b,e,p,q))\\ prop(b,e,p,q):((cond(b,e,p,q)\implies\gcd(p,q)是s_{[b,e)}的周期) \]
总的思路就是将$(b,e,p,q)$拆解为两个子问题$(b,e-p,q-p,p)$和$(b+p,e,q-p,p)$。

首先证明$cond(b,e,p,q)\implies(cond(b,e-\min\{p,q\},|q-p|,\min\{p,q\})\wedge cond(b+\min\{p,q\},e,|q-p|,\min\{p,q\}))$，以便后面将$prop$拆解。不妨设$p\leq q$（如果不成立，交换两者即可）。记$d=q-p$。于是等价于要证

\[cond(b,e,p,q)\implies(cond(b,e-p,d,p)\wedge cond(b+p,e,d,p)) \]
根据更相减损术，$\gcd(p,q)=\gcd(q-p,p)=\gcd(d,p)$，所以

\[\begin{align*} e-(b+p)=(e-p)-b&=(e-b)-p\\ &\geq(p+q)-\gcd(p,q)-p\\ &=((q-p)+p)-\gcd(q-p,p)\\ &=(d+p)-\gcd(d,p) \end{align*} \]
所以$cond_0(b,e-p,d,p)$和$cond_0(b+p,e,d,p)$都得证。

因为$d\geq\gcd(d,p)=\gcd(p,q)$，所以

\[e-b\geq(p+q)-\gcd(p,q)=2p+(r-\gcd(p,q))\geq2p\tag{$*$} \]
因为$p,q$均为$s_{[b,e)}$的周期，所以

\[\forall i\in\mathbb N\cap[b,e-p)\ (s_i=s_{i+p})\\ \forall i\in\mathbb N\cap[b,e-q)\ (s_i=s_{i+q}) \]
所以

\[\forall i\in\mathbb N\cap[b,e-p-d)\ (i+r,i+p\leq i+q<e-p-d+q=e-d+d=e) \]
这保证了以下结论中下标不会越界：

\[\forall i\in\mathbb N\cap[b,e-p-d)\ (s_i=s_{i+q}=s_{i+q-p}=s_{i+d}) \]
即$d$是$s_{[b,e-p)}$的周期。而$p$为周期对于$s_{[b,e-p)}$仍然成立，只要$e-p-b\geq p$即$e-b\geq2p$，而这个不等式在$(*)$处已经得证。所以$cond_1(b,e-p,d,p)$得证。

同理有

\[\forall i\in\mathbb N\cap[b+p,e-d)\ (b\leq i-p,i+d<e) \]
所以以下结论中下标不会越界：

\[\forall i\in\mathbb N\cap[b+p,e-d)\ (s_i=s_{i-p}=s_{i-p+q}=s_{i+d}) \]
即$d$是$s_{[b+p,e)}$的周期。同理，$p$为周期对于$s_{[b+p,e)}$也仍然成立。所以$cond_1(b+p,e,d,p)$得证。

综上，$cond(b,e,p,q)\implies(cond(b,e-p,d,p)\wedge cond(b+p,e,d,p))$。

接下来我们证明$(prop(b,e-\min\{p,q\},|q-p|,\min\{p,q\})\wedge prop(b+\min\{p,q\},e,|q-p|,\min\{p,q\}))\implies prop(b,e,p,q)$。

仍然不妨设$p\leq q$，设$d=q-p$，所以等价于要证

\[(prop(b,e-p,q-p,p)\wedge prop(b+p,e,q-p,p))\implies prop(b,e,p,q) \]

逻辑

这里的逻辑有一点复杂：当前命题$prop:cond\implies conclusion$，两个子命题$prop_i:cond_i\implies conclusion_i\ (i\in\{0,1\})$。我们要证明

\[(prop_0\wedge prop_1)\implies prop \]
即证明

\[\forall i\in\{0,1\}\ (cond\implies cond_i\implies conclusion_i)\\ conclusion_0\wedge conclusion_1\implies conclusion \]
其中第一行就是我们之前证明的结果，接下来我们要证明第二行。

如果两个子命题成立，则$\gcd(p,q)$同时是$s_{[b,e-p)}$和$s_{[b+p,e)}$的周期，也就是

\[\forall i\in\mathbb N\cap([b,e-p-\gcd(p,q))\cup[b+p,e-\gcd(p,q)))\ (s_i=s_{i+g}) \]
因为

\[e-b\geq(p+q)+\gcd(p,q)\geq2p+\gcd(p,q) \]
所以

\[e-p-\gcd(p,q)\geq b+p \]
也就是

\[[b,e-p-\gcd(p,q))\cup[b+p,e-\gcd(p,q))=[b,e-\gcd(p,q)) \]
所以

\[\forall i\in\mathbb N\cap[b,e-\gcd(p,q))\ (s_i=s_{i+g}) \]
即$\gcd(p,q)$是$s_{[b,e)}$的周期，得证。

现在我们已经证明了递归证明的下行和上行部分，只需证明递归边界满足条件即可。由更相减损术可知，递归过程中问题的规模一定是严格递减的，且必将到达递归边界（设为第$x$层）$p_x=0,\ q_x=\gcd(p_0,q_0)$，此时$\gcd(p_x,q_x)=q_x$是$s_{[b_x,e_x)}$的周期，$prop(b_x,e_x,p_x,q_x)$成立。

到此，由初始条件不断递归直至边界，再不断由子命题成立得到上一层的命题成立，我们终于证明了这个定理。不容易啊。

设$\Sigma$为字符集，$p_\mathrm m$是字符串$s\in\Sigma^n\ (n\in\mathbb N)$的最小正周期，且$p_\mathrm m\nmid n$。我们需要证明不存在其他任何$\vec s$的周期$p\not\in\{p_\mathrm m,n\}$满足$p\mid n$。

假设有$\vec s$的周期$p\not\in\{p_\mathrm m,n\}$满足$p\mid n$，由$1<\frac{n}{p}\in\mathbb N^*$知$p\leq n/2$。因为$p_\mathrm m$是最小正周期，所以$p_\mathrm m<p\leq n/2$，于是

\[(p_\mathrm m+p)-\gcd(p_\mathrm m,p)\leq n-\gcd(p_\mathrm m,p)\leq n-1<n \]
符合引理条件，所以$\gcd(p_\mathrm m,p)$也是$\vec s$的周期。还是因为$p_\mathrm m$是最小正周期，所以$\gcd(p_\mathrm m,p)\geq p_\mathrm m$，所以$\gcd(p_\mathrm m,p)=p_\mathrm m$，即$p_\mathrm m\mid p$。因为$p_\mathrm m\nmid n$，所以$p\nmid n$，与假设矛盾，故假设不成立。

字符串哈希

设$\Sigma$为字符集。定义一个字符串$\vec s\in\Sigma^n$的哈希函数为

\[h(\vec s)=\left(\sum_{i=0}^{n-1}s_i\cdot B^{n-1-i}\right)\bmod M \]

其中$B,M$是预先设定的参数，通常都取质数，且满足$\max\Sigma<B<M$。

背诵：$B=233,\ M=993244853$。注意，一个常见的模数是$998244353$，这里只是一个长得比较像的玩意。

可以通过递推法计算$\vec s$所有前缀的哈希值：

\[\forall i\in\mathbb N\cap[0,n)\ (h(s[0,i+1))=(h(s[0,i))\cdot B+s_i)\bmod M) \]

利用所有前缀的哈希值可以快速计算任意子串的哈希值：

\[(\forall i,l\in\mathrm N,\ i+l\leq n)\ (h(s[i,i+l))=(h(s[0,i+l))-h(s[0,i))\cdot B^l)\bmod M) \]

不过需要注意负数取模的问题。

Manacher算法

OI-Wiki的解释相当详细清晰。

/**
 * @return {odd, even}, both denote the number of palindrome subsequences
 * 	centering around each elements (consider the right one as the center for
 * 	even-length palindromes)
 */
template <typename Iter, typename = RequireFwdIter<Iter>>
inline std::pair<std::vector<size_t>, std::vector<size_t>>
manacher(Iter begin, Iter end) {
	size_t n = std::distance(begin, end);
	if (!n) {
		return {{}, {}};
	}
	auto odd = std::vector<size_t>(n), even = std::vector<size_t>(n);
	odd[0] = 1;
	for (size_t i = 1, b = 0, e = 1; i != n; ++i) {
		if (i < e) {
			odd[i] = std::min(odd[b + e - 1 - i], e - i);
		}
		while (odd[i] <= i && i + odd[i] < n &&
			   *std::next(begin, i - odd[i]) ==
				   *std::next(begin, i + odd[i])) {
			++odd[i];
		}
		if (i + odd[i] > e) {
			b = i - odd[i] + 1;
			e = i + odd[i];
		}
	}
	for (size_t i = 0, b = 0, e = 0; i != n; ++i) {
		if (i < e) {
			even[i] = std::min(even[b + e - i], e - i);
		}
		while (even[i] < i && i + even[i] < n &&
			   *std::next(begin, i - even[i] - 1) ==
				   *std::next(begin, i + even[i])) {
			++even[i];
		}
		if (i + even[i] > e) {
			b = i - even[i];
			e = i + even[i];
		}
	}
	return {odd, even};
}

时间复杂度：$O(n)$。
空间复杂度：$O(1)$额外空间。

[2025-08-07] SA、Trie和PAM

后缀数组

设$\Sigma$为字符集，字符串$\vec s\in\Sigma^n\ (n\in\mathbb N)$。称数组$\vec{sa}\in\mathbb N^n$为字符串$\vec s$的后缀数组，当且仅当其满足对任意$i\in\mathbb N\cap[0,n)$，$s[sa_i,n)$是$\vec s$所有后缀中字典序第$i$小的。求后缀数组的同时，记录排名数组$\vec{ra}$，使得对任意$i\in\mathbb N\cap[0,n)$，$ra_i$是$s[i,n)$在$\vec s$所有后缀中按字典序排序后的排名，因此有

\[sa_{ra_i}=i=ra_{sa_i} \]

朴素方法

求后缀数组的朴素方法就是直接枚举所有的后缀然后排序。因为比较两个后缀字典序的时间复杂度为$O(n)$，所以总的时间复杂度为$O(n^2\log n)$，空间复杂度为$O(n^2)$。

第一个优化

利用倍增思想优化。详见OI-Wiki。

template <typename Iter>
inline std::pair<std::vector<size_t>, std::vector<size_t>>
sufArrOf(Iter begin, Iter end) {
	std::vector<size_t> ra, old_ra;
	discretize(begin, end, ra);
	size_t n = ra.size();
	auto sa = std::vector<size_t>(n);

	for (size_t i = 0; i != n; ++i) {
		sa[i] = i;
	}

	for (size_t len = 1; len < n; len <<= 1) {
		old_ra = ra;
		auto cmp = [&](size_t lhs, size_t rhs) {
			if (old_ra[lhs] != old_ra[rhs]) {
				return old_ra[lhs] < old_ra[rhs];
			}
			if (lhs + len >= n) {
				return (rhs + len < n);
			}
			return (rhs + len < n && old_ra[lhs + len] < old_ra[rhs + len]);
		};
		std::sort(sa.begin(), sa.end(), cmp);
		ra[sa[0]] = 0;
		for (size_t i = 1; i != n; ++i) {
			ra[sa[i]] = ra[sa[i - 1]] + cmp(sa[i - 1], sa[i]);
		}
		if (ra[sa.back()] == n - 1) {
			break;
		}
	}

	return {sa, ra};
}

对于字符串，其实可以

for (size_t i = 0; i != s.size(); ++i) {
	sa[i] = i;
	ra[i] = s[i];
}

但是这里为了使其适用于任意定义了全序关系的序列，我使用了discretize。

时间复杂度：$O(n(\log n)^2)$。
空间复杂度：$O(n)$。

第二个优化

我们需要进行双关键字排序，所以可以使用基数排序（radix sort）；而关键字的值域是$\mathbb N\cap[0,n)$，所以可以每一轮排序可以用计数排序（counting sort）实现。

template <typename Iter>
inline std::pair<std::vector<size_t>, std::vector<size_t>>
sufArrOf(Iter begin, Iter end) {
	std::vector<size_t> ra;
	size_t unq = discretize(begin, end, ra).size(), n = ra.size();
	auto cnt = std::vector<size_t>(n);
	auto sa = std::vector<size_t>(n),
		 tmp = std::vector<size_t>(n);

	for (size_t i = 0; i < n; ++i) {
		++cnt[ra[i]];
	}
	for (size_t i = 1; i < unq; ++i) {
		cnt[i] += cnt[i - 1];
	}
	for (size_t i = n; (i--) > 0;) {
		sa[--cnt[ra[i]]] = i;
	}

	for (size_t len = 1, tot; len < n; len <<= 1) {
		// Sort sa so that ra[sa[i] + len] <= ra[sa[i + 1] + len]
		tot = 0;
		for (size_t i = n - len; i != n; ++i) {
			tmp[tot++] = i;
		}
		for (size_t i = 0; i != n; ++i) {
			if (sa[i] >= len) {
				tmp[tot++] = sa[i] - len;
			}
		}

		// Stably sort sa so that ra[sa[i]] <= ra[sa[i + 1]]
		std::fill(cnt.begin(), cnt.begin() + unq, 0);
		for (auto i : tmp) {
			++cnt[ra[i]];
		}
		for (size_t i = 1; i != unq; ++i) {
			cnt[i] += cnt[i - 1];
		}
		for (size_t i = n; (i--) > 0;) {
			sa[--cnt[ra[tmp[i]]]] = tmp[i];
		}

		// Update ra
		std::copy(ra.begin(), ra.end(), tmp.begin());
		ra[sa[0]] = 0;
		tot = 1;
		for (size_t i = 1; i != n; ++i) {
			if (tmp[sa[i]] == tmp[sa[i - 1]] &&
				(sa[i] + len < n) == (sa[i - 1] + len < n) &&
				((sa[i] + len < n)
					 ? (tmp[sa[i] + len] == tmp[sa[i - 1] + len])
					 : true)) {
				ra[sa[i]] = tot - 1;
			} else {
				ra[sa[i]] = (tot++);
			}
		}
		if ((unq = tot) == n) {
			break;
		}
	}

	return {sa, ra};
}

时间复杂度：$O(n\log n)$。
空间复杂度：$O(n)$。

据说有$O(n)$的做法，但是太复杂了，而且一般没必要。

height数组

令$lcp(i,j)$表示字符串$\vec s$分别从$i$和$j$开始的后缀的最长公共前缀（longest common prefix）的长度。height数组$\vec h\in\mathbb N^n$定义为

\[\vec h=(lcp(sa_{i-1},sa_i))_{i=0}^{n-1} \]

即$h_i$为排名第$i$位的后缀与排在它前面的后缀的LCP长度。特别地，规定$h_0=0$。于是有

\[h_{ra_i}=lcp(sa_{ra_i},sa_{ra_i-1})=lcp(i,sa_{ra_i-1}) \]

即$h_{ra_i}$是从$i$处开始的后缀与排在它前一位的后缀的LCP长度。

最长公共前缀引理（LCP lemma）：$h_{ra_i}\geq h_{ra_{i-1}}-1$。

证明

OI-Wiki的证明相当易懂，我就不再抄写一遍了。

计算height数组的Kasai算法

利用LCP引理就可以以$O(n)$的时间复杂度从$\vec{sa},\vec{ra}$导出$\vec h$。以下代码中len维护当前LCP长度。每当i移动到下一位时，先令--len，然后再暴力向后匹配即可。

for (size_t i = 0, j, len = 0; i != n; ++i) {
	if (!ra[i]) {
		len = 0;
		continue;
	}
	if (len) {
		--len;
	}
	for (j = sa[ra[i] - 1];
		 *std::next(begin, i + len) == *std::next(begin, j + len);
		 ++len)
		;
	height[ra[i]] = len;
}

时间复杂度：$O(n)$。这是因为在暴力匹配的过程中，len始终是自增的，而len = 0只会执行一次，这意味着++len的总执行次数不超过$2n$。
空间复杂度：$O(1)$额外空间。

最终代码：

/**
 * @return {sa, ra, height}
 */
template <typename Iter>
inline std::tuple<std::vector<size_t>, std::vector<size_t>, std::vector<size_t>>
sufArrOf(Iter begin, Iter end) {
	std::vector<size_t> ra;
	size_t unq = discretize(begin, end, ra).size(), n = ra.size();
	auto cnt = std::vector<size_t>(n);
	auto sa = std::vector<size_t>(n),
		 height = std::vector<size_t>(n),
		 tmp = std::vector<size_t>(n);

	for (size_t i = 0; i < n; ++i) {
		++cnt[ra[i]];
	}
	for (size_t i = 1; i < unq; ++i) {
		cnt[i] += cnt[i - 1];
	}
	for (size_t i = n; (i--) > 0;) {
		sa[--cnt[ra[i]]] = i;
	}

	for (size_t len = 1, tot; len < n; len <<= 1) {
		// Sort sa so that ra[sa[i] + len] <= ra[sa[i + 1] + len]
		tot = 0;
		for (size_t i = n - len; i != n; ++i) {
			tmp[tot++] = i;
		}
		for (size_t i = 0; i != n; ++i) {
			if (sa[i] >= len) {
				tmp[tot++] = sa[i] - len;
			}
		}

		// Stably sort sa so that ra[sa[i]] <= ra[sa[i + 1]]
		std::fill(cnt.begin(), cnt.begin() + unq, 0);
		for (auto i : tmp) {
			++cnt[ra[i]];
		}
		for (size_t i = 1; i != unq; ++i) {
			cnt[i] += cnt[i - 1];
		}
		for (size_t i = n; (i--) > 0;) {
			sa[--cnt[ra[tmp[i]]]] = tmp[i];
		}

		// Update ra
		std::copy(ra.begin(), ra.end(), tmp.begin());
		ra[sa[0]] = 0;
		tot = 1;
		for (size_t i = 1; i != n; ++i) {
			if (tmp[sa[i]] == tmp[sa[i - 1]] &&
				(sa[i] + len < n) == (sa[i - 1] + len < n) &&
				((sa[i] + len < n)
					 ? (tmp[sa[i] + len] == tmp[sa[i - 1] + len])
					 : true)) {
				ra[sa[i]] = tot - 1;
			} else {
				ra[sa[i]] = (tot++);
			}
		}
		if ((unq = tot) == n) {
			break;
		}
	}

	// Calculate height
	for (size_t i = 0, j, len = 0; i != n; ++i) {
		if (!ra[i]) {
			len = 0;
			continue;
		}
		if (len) {
			--len;
		}
		for (j = sa[ra[i] - 1];
			 *std::next(begin, i + len) == *std::next(begin, j + len);
			 ++len)
			;
		height[ra[i]] = len;
	}

	return {sa, ra, height};
}

[2025-08-10] ACAM和SAM

[2025-08-10] 强连通分量、双连通分量、割点和割边

Tarjan算法

Tarjan算法适合处理连通性相关的问题，用途相当广泛。它的核心是在图$G=(V,E)$上DFS的过程中维护两个数组$\overrightarrow{dfn}\in\mathbb N^{|V|}$（depth-first number）和$\overrightarrow{low}\in\mathbb N^{|V|}$（low-link value），对于$i\in\mathbb N\cap[0,|V|)$，$dfn_i$是DFS过程中第一次访问$i$的时间（或者说，顺序），$low_i$则表示从顶点$i$出发沿未走过的边（有向图是所有从节点$i$引出的边，无向图则是所有与$i$相连的边中除从父节点走到$i$的那一条以外的所有边）能够到达的顶点中（包括它自己）最小的$dfn$。

\[\begin{array}{rl} 1 & \textbf{function }tarjan(V,E)\\ 2 & \qquad tm\gets0\\ 3 & \qquad\overrightarrow{dfn}\gets(-1)_{i=0}^{|V|-1}\\ 4 & \qquad\textbf{function }dfs(pre,cur)\\ 5 & \qquad\qquad dfn_{cur}\gets tm\\ 6 & \qquad\qquad low_{cur}\gets dfn_{cur}\\ 7 & \qquad\qquad tm\gets tm+1\\ 8 & \qquad\qquad\textbf{for each }(u, v)\textbf{ in } E \\ 9 & \qquad\qquad\qquad nxt\gets\begin{cases}v&(u=cur)\\u&(u\neq cur)\end{cases}\\ 10 & \qquad\qquad\textbf{if }dfn_{nxt}=-1\\ 11 & \qquad\qquad\qquad dfs(cur,nxt)\\ 12 & \qquad\qquad\qquad low_{cur}\gets\min\{low_{cur},low_{nxt}\}\\ 13 & \qquad\qquad\textbf{else if }nxt\neq pre\\ 14 & \qquad\qquad\qquad low_{cur}\gets\min\{low_{cur},dfn_{nxt}\}\\ 15 & \qquad\textbf{for each }i\textbf{ in }\mathbb N\cap[0,|V|)\\ 16 & \qquad\qquad\textbf{if }dfn_i=-1\\ 17 & \qquad\qquad\qquad dfs(i,i)\\ 18 & \qquad\textbf{return }(\overrightarrow{dfn},\overrightarrow{low}) \end{array} \]

下面的Tarjan算法都是在此DFS的基础上增添内容。

强连通分量

对于有向图$G=(V,E)$，我们称其强连通，当且仅当对任意$u,v\in V$，$u$与$v$双向可达。

称一个有向图的极大连通子图为其强连通分量（Strongly Connected Components, SCC）。

	std::vector<std::vector<size_t>> tarjanSccs() const {
		static_assert(is_directed,
					  "Tarjan's algorithm for strongly connected components "
					  "is only applicable to directed graphs.");
		std::vector<size_t> stk;
		stk.reserve(m_adj.size());
		std::vector<bool> in_stk(m_adj.size());
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> sccs;
		size_t tm = 0;
		std::function<void(size_t)> dfs = [&](size_t cur) -> void {
			stk.emplace_back(cur);
			in_stk[cur] = true;
			low[cur] = (dfn[cur] = (tm++));
			for (auto i : m_adj[cur]) {
				auto &edge = m_edges[i];
				if (edge.u != cur) { // edges pointing to cur
					continue;
				}
				if (dfn[edge.v] == size_t(-1)) {
					dfs(edge.v);
					minEq(low[cur], low[edge.v]);
				} else if (in_stk[edge.v]) {
					minEq(low[cur], dfn[edge.v]);
				}
			}
			if (dfn[cur] == low[cur]) {
				sccs.emplace_back();
				size_t vert;
				do {
					sccs.back().emplace_back(vert = stk.back());
					stk.pop_back();
					in_stk[vert] = false;
				} while (vert != cur);
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				dfs(i);
			}
		}
		return sccs;
	}

时间复杂度：$O(|V|+|E|)$。
空间复杂度：$O(|V|)$。

割点

对于连通无向图$G=(V,E)$，若存在$v\in V$满足$(V\backslash\{v\}, E)$不连通，则称$v$为$G$的一个割点（cut vertex）。

仍然在DFS生成树上考虑，一个点$cur$是割点当且仅当存在它的一个子节点$nxt$，使得$low_{nxt}\geq dfn_{cur}$，即$nxt$仅通过$cur$这个点与其他点相连。这里取等是可以的，因为即便有很多条重边$(cur,nxt)$，都不影响以$nxt$为根的子树对于$cur$的依赖。

特别地，当$cur$是DFS生成树的根节点时，$cur$是割点当且仅当$cur$有超过一个子树。

证明

充分性：当$cur$有超过一个子树时，它们都依赖$cur$来相互沟通。如果它们之间存在别的沟通方式，就不会形成不同的子树，而是会在同一次向下搜索的过程中都访问到。

必要性：如果$cur$没有子树，那显然不是割点。如果有且仅有一颗子树，那去掉$cur$对于其他点的连通性也没有影响。

之所以要特判，是因为对于根节点来说，必有$low_{nxt}\geq dfn_{cur}$，而上面的证明说明，此时$cur$不一定是割点。

以下函数求的是一个无向图中各个强连通分量的所有割点。

std::vector<std::vector<size_t>> tarjanCutVerticesOfSccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for vertex cut "
					  "is only applicable to undirected graphs.");
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> cut_verts;
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) -> void {
			low[cur] = (dfn[cur] = (tm++));
			size_t children = 0;
			bool not_in_cut = true;
			for (auto i : m_adj[cur]) {
				size_t nxt = ((m_edges[i].u == cur) ? m_edges[i].v : m_edges[i].u);
				if (dfn[nxt] == size_t(-1)) {
					++children;
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (not_in_cut && ~pre && low[nxt] >= dfn[cur]) {
						not_in_cut = false;
						cut_verts.back().emplace_back(cur);
					}
				} else if (nxt != pre) {
					minEq(low[cur], dfn[nxt]);
				}
			}
			if (not_in_cut && pre == size_t(-1) && children > 1) {
				cut_verts.back().emplace_back(cur);
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				cut_verts.emplace_back();
				dfs(-1, i);
			}
		}
		return cut_verts;
	}

时间复杂度：$O(|V|+|E|)$。
空间复杂度：$O(|V|)$。

割边

对于连通无向图$G=(V,E)$，若存在$e\in E$满足$(V,E\backslash\{e\})$不连通，则称$v$为$G$的一个割边/桥（bridge）。

类似于上面求割点的想法，在DFS生成树上，一个节点$cur$与它的一个子节点$nxt$之间的边是割边，当且仅当$low_{nxt}>dfn_{cur}$，即$nxt$没有返祖边。这时以$nxt$为根的子树与其他部分连通的唯一途径就是边$(cur,nxt)$。

不过如果有重边的话，这里的判断就会出问题，因为我们之前都是通过判断$nxt$是否等于$pre$来决定$(cur,nxt)$是不是走过的边，而当$cur$与$nxt$之间有重边时，这个判断就是错误的，重边会被忽略掉，从而使本不是割边的$(cur,nxt)$（因为割掉一条，还有重边）也被误判为割边了。一个简单的方案是只允许$nxt=pre$一次（下面用has_edge_to_pre记录），下一次等于的时候就要判断为返祖边了。

以下函数求的是一个无向图中各个强连通分量的所有割边。has_edge_to_pre用于处理有重边的情况。

	std::vector<std::vector<size_t>> tarjanBridgesOfSccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for bridges "
					  "is only applicable to undirected graphs.");
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> bridges;
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) -> void {
			low[cur] = (dfn[cur] = (tm++));
			bool has_edge_to_pre = false;
			for (auto i : m_adj[cur]) {
				size_t nxt = ((m_edges[i].u == cur) ? m_edges[i].v : m_edges[i].u);
				if (dfn[nxt] == size_t(-1)) {
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (low[nxt] > dfn[cur]) {
						bridges.back().emplace_back(i);
					}
				} else if (nxt != pre || has_edge_to_pre) {
					minEq(low[cur], dfn[nxt]);
				} else {
					has_edge_to_pre = true;
				}
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				bridges.emplace_back();
				dfs(-1, i);
			}
		}
		return bridges;
	}

时间复杂度：$O(|V|+|E|)$。
空间复杂度：$O(|V|)$。

可以将求割点和割边的部分整合到一起，提高代码复用率：

	std::pair<std::vector<std::vector<size_t>>, std::vector<std::vector<size_t>>>
	tarjanCutVerticesAndBridgesOfSccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for cut vertices and bridges "
					  "is only applicable to undirected graphs.");
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> cut_verts, bridges;
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) -> void {
			low[cur] = (dfn[cur] = (tm++));
			size_t children = 0;
			bool not_in_cut = true, has_edge_to_pre = false;
			for (auto i : m_adj[cur]) {
				size_t nxt = ((m_edges[i].u == cur) ? m_edges[i].v : m_edges[i].u);
				if (dfn[nxt] == size_t(-1)) {
					++children;
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (not_in_cut && ~pre && low[nxt] >= dfn[cur]) {
						not_in_cut = false;
						cut_verts.back().emplace_back(cur);
					}
					if (low[nxt] > dfn[cur]) {
						bridges.back().emplace_back(i);
					}
				} else if (nxt != pre || has_edge_to_pre) {
					minEq(low[cur], dfn[nxt]);
				} else {
					has_edge_to_pre = true;
				}
			}
			if (not_in_cut && pre == size_t(-1) && children > 1) {
				cut_verts.back().emplace_back(cur);
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				cut_verts.emplace_back();
				bridges.emplace_back();
				dfs(-1, i);
			}
		}
		return {cut_verts, bridges};
	}

点双连通分量

在无向图$G=(V,E)$中，若点$v_0,v_1$满足$\forall v\in V\backslash\{v_0,v_1\}\ (v_0与v_1在图(V\backslash\{v\},E)中连通)$，则称$v_0$与$v_1$点双连通，若任意两点均点双连通，则称$G$点双连通（2-vertex-/biconnected）。

对于无向图$G$，称其一个极大的点双连通子图为$G$的一个点双连通分量（2-vertex/biconnected component/block）。

因为一个点双连通分量内部是没有割点的，所以点双连通分量之间要么完全不连通，要么就一定有共同的割点，所以可以在求割点的Tarjan算法基础上，将割点的各个子树分别并上这个割点之后划分为双连通分量。具体地，在DFS生成树上，如果点$cur$是割点，$nxt$是$cur$的一个子节点，则以$nxt$为根的子树并上$\{cur\}$就是原图的一个点双连通分量。仿照求强连通分量的过程，我们用一个栈来记录当前访问过且可以被加入点双连通分量的顶点。根节点同样需要特殊处理，且需要在主循环中将栈中剩余的点都划入一个新的点双连通分量，细节很多。

	std::vector<std::vector<size_t>> tarjanVbccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for 2-vertex-connected components "
					  "is only applicable to undirected graphs.");
		std::vector<size_t> stk;
		stk.reserve(m_adj.size());
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> vbccs;
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) -> void {
			dfn[cur] = low[cur] = tm++;
			stk.push_back(cur);
			size_t children = 0;
			for (auto i : m_adj[cur]) {
				size_t nxt = (m_edges[i].u == cur)
								 ? m_edges[i].v
								 : m_edges[i].u;
				if (dfn[nxt] == size_t(-1)) {
					++children;
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (~pre ? (low[nxt] >= dfn[cur]) : (children > 1)) {
						vbccs.emplace_back(1, cur);
						size_t vert;
						do {
							vbccs.back().push_back(vert = stk.back());
							stk.pop_back();
						} while (vert != nxt);
					}
				} else if (nxt != pre) {
					minEq(low[cur], dfn[nxt]);
				}
			}
		};
		for (size_t i = 0; i < m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				dfs(-1, i);
				if (!stk.empty()) {
					vbccs.emplace_back(std::move(stk));
				}
			}
		}
		return vbccs;
	}

时间复杂度：$O(|V|+|E|)$。
空间复杂度：$O(|V|)$。

边双连通分量

在无向图$G=(V,E)$中，若点$v_0,v_1$满足$\forall e\in E\ (v_0与v_1在图(V,E\backslash\{e\})中连通)$，则称$v_0$与$v_1$点双连通，若任意两点均点双连通，则称$G$边双连通（2-edge-connected）。

对于无向图$G$，称其一个极大的边双连通子图为$G$的一个边双连通分量（2-edge-connected component）。

类似于求点双连通分量，我们只需求出割边，然后将割边远离DFS生成树的树根的那一端划入边双连通分量即可。

	std::vector<std::vector<size_t>> tarjanEbccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for 2-edge-connected components "
					  "is only applicable to undirected graphs.");
		std::vector<size_t> stk;
		stk.reserve(m_adj.size());
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		std::vector<std::vector<size_t>> ebccs;
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) -> void {
			dfn[cur] = low[cur] = tm++;
			stk.emplace_back(cur);
			bool has_edge_to_pre = false;
			for (auto i : m_adj[cur]) {
				size_t nxt = (m_edges[i].u == cur)
								 ? m_edges[i].v
								 : m_edges[i].u;
				if (dfn[nxt] == size_t(-1)) {
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (low[nxt] > dfn[cur]) {
						ebccs.emplace_back();
						size_t vert;
						do {
							ebccs.back().emplace_back(vert = stk.back());
							stk.pop_back();
						} while (vert != nxt);
					}
				} else if (nxt != pre || has_edge_to_pre) {
					minEq(low[cur], dfn[nxt]);
				} else {
					has_edge_to_pre = true;
				}
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				dfs(-1, i);
				if (stk.size()) {
					ebccs.emplace_back(std::move(stk));
				}
			}
		}
		return ebccs;
	}

时间复杂度：$O(|V|+|E|)$。

空间复杂度：$O(|V|)$。

因为要求双连通分量，就要求割点或割边，所以完全可以将它们全都融合起来，狠狠复用代码。

	/**
	 * @return {{cut_verts, vbccs}, {bridges, ebccs}}
	 */
	inline std::pair<std::pair<std::vector<std::vector<size_t>>,
							   std::vector<std::vector<size_t>>>,
					 std::pair<std::vector<std::vector<size_t>>,
							   std::vector<std::vector<size_t>>>>
	tarjanCutAndBccs() const {
		static_assert(!is_directed,
					  "Tarjan's algorithm for cut vertices, bridges "
					  "and biconnected components is only applicable to "
					  "undirected graphs.");
		std::vector<std::vector<size_t>> cut_verts, vbccs, bridges, ebccs;
		std::vector<size_t> vbcc_stk, ebcc_stk;
		std::vector<size_t> dfn(m_adj.size(), size_t(-1)), low(m_adj.size());
		size_t tm = 0;
		std::function<void(size_t, size_t)> dfs = [&](size_t pre, size_t cur) {
			low[cur] = (dfn[cur] = (tm++));
			vbcc_stk.emplace_back(cur);
			ebcc_stk.emplace_back(cur);
			size_t children = 0;
			bool not_in_cut = true, has_edge_to_pre = false;
			for (auto i : m_adj[cur]) {
				size_t nxt = ((m_edges[i].u == cur)
								  ? m_edges[i].v
								  : m_edges[i].u);
				if (dfn[nxt] == size_t(-1)) {
					++children;
					dfs(cur, nxt);
					minEq(low[cur], low[nxt]);
					if (~pre ? (low[nxt] >= dfn[cur]) : (children > 1)) {
						if (not_in_cut) {
							not_in_cut = false;
							cut_verts.back().emplace_back(cur);
						}
						vbccs.emplace_back(1, cur);
						size_t vert;
						do {
							vbccs.back().push_back(vert = vbcc_stk.back());
							vbcc_stk.pop_back();
						} while (vert != nxt);
					}
					if (low[nxt] > dfn[cur]) {
						bridges.back().emplace_back(i);
						ebccs.emplace_back();
						size_t vert;
						do {
							ebccs.back().emplace_back(vert = ebcc_stk.back());
							ebcc_stk.pop_back();
						} while (vert != nxt);
					}
				} else if (nxt != pre || has_edge_to_pre) {
					minEq(low[cur], dfn[nxt]);
				} else {
					has_edge_to_pre = true;
				}
			}
		};
		for (size_t i = 0; i != m_adj.size(); ++i) {
			if (dfn[i] == size_t(-1)) {
				cut_verts.emplace_back();
				bridges.emplace_back();
				dfs(-1, i);
				if (vbcc_stk.size()) {
					vbccs.emplace_back(std::move(vbcc_stk));
				}
				if (ebcc_stk.size()) {
					ebccs.emplace_back(std::move(ebcc_stk));
				}
			}
		}
		return {{cut_verts, vbccs}, {bridges, ebccs}};
	}

爽！

posted @ 2025-08-16 03:07 我就是蓬蒿人阅读(44) 评论(0) 收藏举报

刷新页面返回顶部

仰天大笑出门去

BITACM 2025 暑期集训 B组 笔记

Notes

[2025-07-15] DP基础

0/1背包

完全背包

多重背包

优化

最长上升子序列（LIS）

[2025-07-16] 树状数组和线段树

树状数组

线段树

[2025-07-17] DFS、BFS、最短路和拓扑排序

最短路

拓扑排序

[2025-07-19] 单调栈和单调队列

单调栈

单调队列

[2025-07-20] 贪心、二分和双指针

贪心

二分

双指针

[2025-07-23] 位运算和状压DP

位运算

状压DP

最短哈密顿回路

[2025-07-26] 稀疏表、左偏树和LCA

稀疏表

[2025-07-31] FFT和NTT

FFT

问题

朴素方法

第一个优化：点表示

第二个优化：函数的奇偶性

拓展到复数域

FFT

IFFT

递归实现FFT/IFFT

位逆序置换优化

蝶形运算优化

非递归实现FFT

NTT

费马小定理

[2025-08-02] 数论基础

乘法逆元

扩展欧几里得算法

快速幂

递推法

质数筛

埃氏筛（Sieve of Eratosthenes）

欧氏/线性筛（Euler’s/Linear Sieve）

[2025-08-03] 狄利克雷卷积、莫反和杜教筛

数论基础

同余类和剩余系

数论函数

积性函数

狄利克雷卷积

莫比乌斯反演

欧拉函数

[2025-08-04] BITACM SC 2025 Rating #3

组合数学

错排公式

排列数和组合数

快速取模

[2025-08-06] KMP、字符串哈希和Manacher

前缀函数

前缀函数的计算

第一个优化

第二个优化

KMP算法

计算序列的最小正周期

字符串哈希

Manacher算法

[2025-08-07] SA、Trie和PAM

后缀数组

朴素方法

第一个优化

第二个优化

height数组

计算height数组的Kasai算法

BITACM 2025 暑期集训 B组笔记