Chuoer47's Blog

Appearance

算法模板 - 数学 ​

本文整理数学相关的算法模板，按难度递增排列：质数 → 快速幂 → 组合数学 → 排列组合 → 线性根 → 高级技巧 → 几何计算。

1. 质数筛 ​

1.1 埃拉托斯特尼筛法（埃筛） ​

核心思想：从 2 开始，把每个质数的倍数标记为合数。剩下的就是质数。

python
N = 1000010
is_prime = [True] * N
primes = []

def eratosthenes():
    for i in range(2, N):
        if is_prime[i]:
            primes.append(i)
            for j in range(i * i, N, i):  # 从 i² 开始标记
                is_prime[j] = False

eratosthenes()
print(primes[:10])  # [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

为什么从 i² 开始？ 小于 i² 的倍数（如 2i, 3i, ...）已经被更小的质数标记过了。比如 6 = 2×3 在 i=2 时就被标记了，不需要在 i=3 时再标记一次。

时间复杂度：O(n log log n)，实际很快。

1.2 欧拉筛（线性筛） ​

核心思想：每个合数只被它的最小质因数标记一次，保证 O(n)。

python
def euler_sieve(n: int) -> list:
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    primes = []

    for i in range(2, n + 1):
        if is_prime[i]:
            primes.append(i)
        for p in primes:
            if p * i > n:
                break
            is_prime[p * i] = False
            if i % p == 0:  # 关键：p 是 i 的最小质因数时停止
                break

    return primes

print(euler_sieve(30))  # [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

为什么 i % p == 0 时停止？ 此时 p 是 i 的最小质因数。如果继续用更大的质数 q 去标记 qi，那么 qi 的最小质因数应该是 p（因为 p < q 且 p | i），应该在后续由 p 去标记。这样保证每个合数只被标记一次。

2. 质因数分解 ​

python
def prime_factors(n):
    """返回 (factors, powers)，如 12 → ([2, 3], [2, 1]) 表示 2² × 3¹"""
    factor, power = [], []

    cnt = 0
    while n % 2 == 0:
        cnt += 1; n //= 2
    if cnt:
        factor.append(2); power.append(cnt)

    for i in range(3, int(n ** 0.5) + 1, 2):
        cnt = 0
        while n % i == 0:
            cnt += 1; n //= i
        if cnt:
            factor.append(i); power.append(cnt)

    if n > 2:
        factor.append(n); power.append(1)

    return factor, power

print(prime_factors(12))  # ([2, 3], [2, 1])

为什么只循环到 √n？ 如果 n 有大于 √n 的因子，那它最多只有一个（否则乘积 > n），就是最后剩余的那个 n。

3. 快速幂 ​

核心思想：把指数拆成二进制，比如 a^13 = a^8 × a^4 × a^1，只需要 log(n) 次乘法。

python
def power(base, exp, mod):
    """计算 base^exp mod mod，O(log exp)"""
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = result * base % mod
        base = base * base % mod
        exp >>= 1
    return result

Python 内置：pow(base, exp, mod) 直接支持模幂运算，底层用 C 实现，比手写快。

4. 矩阵快速幂 ​

适用场景：线性递推的第 n 项（如斐波那契数列），用矩阵快速幂可以 O(k³ log n) 求解。

核心思想：把递推关系写成矩阵形式 [F(n), F(n-1)] = [F(1), F(0)] × M^(n-1)，然后用快速幂求 M^(n-1)。

python
class Matrix:
    def __init__(self, data: list):
        self.data = data
        self.rows = len(data)
        self.cols = len(data[0]) if self.rows else 0

    def __mul__(self, other):
        result = [[0] * other.cols for _ in range(self.rows)]
        for i in range(self.rows):
            for j in range(other.cols):
                for k in range(self.cols):
                    result[i][j] += self.data[i][k] * other.data[k][j]
        return Matrix(result)

    def power(self, n: int):
        identity = [[1 if i == j else 0 for j in range(self.rows)]
                     for i in range(self.rows)]
        result = Matrix(identity)
        current = self
        while n > 0:
            if n & 1:
                result = result * current
            current = current * current
            n >>= 1
        return result

# 斐波那契：[F(n), F(n-1)] = [F(1), F(0)] × [[1,1],[1,0]]^(n-1)
fib_matrix = Matrix([[1, 1], [1, 0]])
result = fib_matrix.power(10)
print(result.data[0][0])  # F(11) = 89

5. 组合数学 ​

5.1 阶乘与逆元预处理 ​

在质数模数下，组合数 C(n, m) = n! / (m! × (n-m)!) 需要用模逆元。

python
MOD = 10 ** 9 + 7
MAXN = 200000

fac = [1] * (MAXN + 1)
inv_f = [1] * (MAXN + 1)

for i in range(1, MAXN + 1):
    fac[i] = fac[i - 1] * i % MOD

inv_f[MAXN] = pow(fac[MAXN], MOD - 2, MOD)  # 费马小定理
for i in range(MAXN, 0, -1):
    inv_f[i - 1] = inv_f[i] * i % MOD

def comb(n, m):
    if m < 0 or n < m:
        return 0
    return fac[n] * inv_f[m] % MOD * inv_f[n - m] % MOD

费马小定理：当 p 是质数时，a^(p-1) ≡ 1 (mod p)，所以 a^(-1) ≡ a^(p-2) (mod p)。Python 的 pow(a, p-2, p) 直接算。

5.2 C(n, m) mod 10（Lucas + CRT） ​

当模数不是质数（如 10 = 2 × 5），需要分别计算 mod 2 和 mod 5，再用中国剩余定理合并。

python
def comb_mod10(n: int, m: int) -> int:
    if m < 0 or m > n: return 0
    if m == 0: return 1

    def mod2(n, m):
        """C(n,m) mod 2：m 的二进制是 n 的子集时为 1"""
        return 1 if (n & m) == m else 0

    def mod5(n, m):
        """C(n,m) mod 5：Lucas 定理，逐位计算"""
        table = [[1,0,0,0,0],[1,1,0,0,0],[1,2,1,0,0],[1,3,3,1,0],[1,4,1,4,1]]
        result = 1
        while n > 0 or m > 0:
            ni, mi = n % 5, m % 5
            if mi > ni: return 0
            result = result * table[ni][mi] % 5
            n //= 5; m //= 5
        return result

    a, b = mod2(n, m), mod5(n, m)
    return (5 * ((a - b) % 2) + b) % 10  # CRT 合并

6. 排列与组合（回溯） ​

6.1 全排列（无重复元素） ​

python
def permute(nums):
    result, used = [], [False] * len(nums)
    def backtrack(path):
        if len(path) == len(nums):
            result.append(path[:]); return
        for i in range(len(nums)):
            if used[i]: continue
            used[i] = True; path.append(nums[i])
            backtrack(path)
            path.pop(); used[i] = False
    backtrack([])
    return result

6.2 多重集排列（含重复元素去重） ​

python
def permute_unique(nums):
    nums.sort()
    result, used = [], [False] * len(nums)
    def backtrack(path):
        if len(path) == len(nums):
            result.append(path[:]); return
        for i in range(len(nums)):
            if used[i]: continue
            # 关键去重：相同元素必须按顺序使用
            if i > 0 and nums[i] == nums[i-1] and not used[i-1]:
                continue
            used[i] = True; path.append(nums[i])
            backtrack(path)
            path.pop(); used[i] = False
    backtrack([])
    return result

print(permute_unique([1, 1, 2]))  # [[1,1,2], [1,2,1], [2,1,1]]

去重条件 not used[i-1] 的含义：如果前一个相同元素还没被使用，说明当前元素"跳过了前一个"，会产生重复排列。

6.3 组合 C(n, k) ​

python
def combine(n, k):
    result = []
    def backtrack(start, path):
        if len(path) == k:
            result.append(path[:]); return
        for i in range(start, n + 1):
            path.append(i)
            backtrack(i + 1, path)
            path.pop()
    backtrack(1, [])
    return result

6.4 组合总和（元素不可重复使用） ​

python
def combination_sum_unique(candidates, target):
    candidates.sort()
    result = []
    def backtrack(start, path, remaining):
        if remaining == 0:
            result.append(path[:]); return
        for i in range(start, len(candidates)):
            if i > start and candidates[i] == candidates[i-1]:
                continue  # 去重
            if candidates[i] > remaining:
                break
            path.append(candidates[i])
            backtrack(i + 1, path, remaining - candidates[i])
            path.pop()
    backtrack(0, [], target)
    return result

7. 线性基（异或） ​

适用场景：求一组数的最大异或和、判断某个数能否被异或表示。

核心思想：维护一组"基"，每个基的最高位互不相同。插入时从高位到低位贪心消去。

python
class XorBasis:
    def __init__(self, n: int):
        self.b = [0] * n

    def insert(self, x: int) -> None:
        b = self.b
        while x:
            i = x.bit_length() - 1
            if b[i] == 0:
                b[i] = x; return
            x ^= b[i]

    def max_xor(self) -> int:
        res = 0
        for i in range(len(self.b) - 1, -1, -1):
            if self.b[i] and (res ^ self.b[i]) > res:
                res ^= self.b[i]
        return res

basis = XorBasis(32)
for x in [3, 5, 7, 9]:
    basis.insert(x)
print(basis.max_xor())  # 15

为什么贪心正确？ 每个基的最高位互不相同。从高位到低位，如果异或当前基能让结果变大，就异或——因为高位的 1 比低位的所有 1 加起来都大。

8. 幺半群滑动窗口 ​

适用场景：滑动窗口内做可结合、有单位元的运算（如 min、max、gcd、sum），O(1) 追加/弹出/查询。

python
from collections import deque

class MonoidWindow:
    def __init__(self, op, e):
        self.op = op; self.e = e
        self.front = deque()
        self.back = deque()

    def _push_back(self, val):
        agg = self.op(self.back[-1][1] if self.back else self.e, val)
        self.back.append((val, agg))

    def _push_front(self, val):
        agg = self.op(val, self.front[0][1] if self.front else self.e)
        self.front.appendleft((val, agg))

    def append(self, val):
        self._push_back(val)

    def pop(self):
        if not self.front:
            while self.back:
                val, _ = self.back.pop()
                self._push_front(val)
        self.front.popleft()

    def query(self):
        if not self.front:
            return self.back[-1][1] if self.back else self.e
        if not self.back:
            return self.front[-1][1]
        return self.op(self.front[-1][1], self.back[-1][1])

# 示例：滑动窗口最小值
w = MonoidWindow(op=min, e=float('inf'))
for x in [3, 1, 4, 1, 5, 9]:
    w.append(x)
print(w.query())  # 1

9. 根号分解（分块） ​

适用场景：当操作按"步长 k"进行时，对 k 的大小分治——k 小的用差分数组批量处理，k 大的直接暴力。

python
from collections import defaultdict
from math import isqrt

def range_multiply_xor(nums, queries):
    """根号分解：区间乘法查询后的异或值"""
    n = len(nums)
    mod = 10**9 + 7
    B = isqrt(len(queries))
    diff = defaultdict(list)

    for l, r, k, v in queries:
        if k < B:
            if not diff[k]:
                diff[k] = [1] * (n + k)
            diff[k][l] = diff[k][l] * v % mod
            r = r - (r - l) % k + k
            diff[k][r] = diff[k][r] * pow(v, -1, mod) % mod
        else:
            for i in range(l, r + 1, k):
                nums[i] = nums[i] * v % mod

    for k, d in diff.items():
        for start in range(k):
            mul = 1
            for i in range(start, n, k):
                mul = mul * d[i] % mod
                nums[i] = nums[i] * mul % mod

    return 0  # 简化示例

根号分解的精髓：√n 是"暴力"和"预处理"的分界线。步长 < √n 的操作种类少，可以按步长预处理；步长 ≥ √n 的操作每次涉及的元素少（≤ √n 个），直接暴力。

10. 几何 — Point 类 ​

2D 点/向量的基础运算，是所有几何算法的基石。

python
from math import hypot, atan2, sin, cos, acos

EPS = 1e-9

class Point:
    __slots__ = ("x", "y")
    def __init__(self, x=0.0, y=0.0):
        self.x = float(x); self.y = float(y)

    def __add__(self, o): return Point(self.x + o.x, self.y + o.y)
    def __sub__(self, o): return Point(self.x - o.x, self.y - o.y)
    def __mul__(self, k): return Point(self.x * k, self.y * k)

    def dot(self, o):   return self.x * o.x + self.y * o.y       # 点积
    def cross(self, o): return self.x * o.y - self.y * o.x       # 叉积
    def length(self):   return hypot(self.x, self.y)
    def dist(self, o):  return hypot(self.x - o.x, self.y - o.y)

    def rotate(self, rad):
        """逆时针旋转 rad 弧度"""
        return Point(self.x*cos(rad) - self.y*sin(rad),
                     self.x*sin(rad) + self.y*cos(rad))

叉积的几何意义：a × b > 0 → b 在 a 的逆时针方向（左转）；< 0 → 顺时针（右转）；= 0 → 共线。

11. 几何 — 直线类 ​

算法专用直线类：将斜率和截距标准化为最简分数，避免浮点精度问题，支持哈希和相等判断。

python
from math import gcd

class Line:
    def __init__(self, p1, p2):
        x1, y1 = p1; x2, y2 = p2
        dx, dy = x2 - x1, y2 - y1
        if dx == 0:       self.slope = (0, 1)       # 垂直线
        elif dy == 0:     self.slope = (1, 0)       # 水平线
        else:
            g = gcd(abs(dy), abs(dx))
            s = (dy // g, dx // g)
            if s[1] < 0: s = (-s[0], -s[1])
            self.slope = s
        # 截距标准化类似，省略...

    def __eq__(self, other):
        return self.slope == other.slope and self.intercept == other.intercept
    def __hash__(self):
        return hash((self.slope, self.intercept))

为什么要标准化？ 直线 y = 2/4 x + 1 和 y = 1/2 x + 1 是同一条直线，但如果不标准化，它们的哈希值不同。

12. 几何 — 线段类 ​

python
class Segment:
    def __init__(self, p1: Point, p2: Point):
        self.p1, self.p2 = p1, p2

    def distance_to_point(self, p: Point) -> float:
        """点到线段的最短距离"""
        vx = self.p2.x - self.p1.x
        vy = self.p2.y - self.p1.y
        wx = p.x - self.p1.x
        wy = p.y - self.p1.y
        c1 = vx * wx + vy * wy
        if c1 <= 0: return p.dist(self.p1)       # 最近点是 p1
        c2 = vx * vx + vy * vy
        if c2 <= c1: return p.dist(self.p2)       # 最近点是 p2
        b = c1 / c2                                # 投影点
        proj = Point(self.p1.x + b*vx, self.p1.y + b*vy)
        return p.dist(proj)

13. 几何 — 圆类 ​

python
import math

class Circle:
    def __init__(self, center: Point, r: float):
        self.center = center; self.r = r

    def contains_point(self, p: Point) -> bool:
        return self.center.dist(p) <= self.r + 1e-12

    def intersects_circle(self, other: "Circle") -> bool:
        d = self.center.dist(other.center)
        return abs(self.r - other.r) <= d <= self.r + other.r

    def intersection_points(self, other: "Circle") -> list:
        """求两个圆的交点（0/1/2 个）"""
        c1, c2 = self.center, other.center
        d = c1.dist(c2)
        if d > self.r + other.r + 1e-12 or d < abs(self.r - other.r) - 1e-12:
            return []
        a = (self.r**2 - other.r**2 + d*d) / (2*d)
        h = math.sqrt(max(0, self.r**2 - a**a))
        x0 = c1.x + a*(c2.x - c1.x)/d
        y0 = c1.y + a*(c2.y - c1.y)/d
        if h == 0: return [Point(x0, y0)]
        rx = -(c2.y - c1.y) * (h/d)
        ry = (c2.x - c1.x) * (h/d)
        return [Point(x0+rx, y0+ry), Point(x0-rx, y0-ry)]

14. 几何 — 凸包（Andrew 算法） ​

核心思想：按 x 坐标排序，分别构建下凸包和上凸包。用叉积判断是否"左拐"，如果不是就弹出栈顶。

python
def convex_hull(points):
    points = sorted(set(points))
    if len(points) <= 2: return points

    def cross(O, A, B):
        return (A[0]-O[0])*(B[1]-O[1]) - (A[1]-O[1])*(B[0]-O[0])

    lower = []
    for p in points:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)

    upper = []
    for p in reversed(points):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)

    return lower[:-1] + upper[:-1]

15. 几何 — 扫描线（矩形面积并） ​

核心思想：将矩形的左右边界作为"事件"，按 x 坐标排序后逐个处理。用一个数组记录每个 y 区间被覆盖的次数。

python
def rectangle_area_union(rectangles):
    """LeetCode 850. 矩形面积 II"""
    MOD = 10**9 + 7
    hbound = sorted({y for rect in rectangles for y in (rect[1], rect[3])})
    m = len(hbound)
    seg = [0] * (m - 1)

    sweep = []
    for i, (x1, y1, x2, y2) in enumerate(rectangles):
        sweep.append((x1, i, 1))   # 左边界进入
        sweep.append((x2, i, -1))  # 右边界离开
    sweep.sort()

    ans = 0; i = 0
    while i < len(sweep):
        j = i
        while j+1 < len(sweep) and sweep[j+1][0] == sweep[i][0]:
            j += 1
        if j+1 == len(sweep): break

        for k in range(i, j+1):
            _, idx, diff = sweep[k]
            for x in range(m-1):
                if rectangles[idx][1] <= hbound[x] and hbound[x+1] <= rectangles[idx][3]:
                    seg[x] += diff

        cover = sum(hbound[k+1] - hbound[k] for k in range(m-1) if seg[k] > 0)
        ans += cover * (sweep[j+1][0] - sweep[j][0])
        i = j + 1

    return ans % MOD

16. Berlekamp-Massey 算法 ​

适用场景：给定数列的前若干项，求最短线性递推式，然后用 Kitamasa 算法 O(k² log n) 求第 n 项。

python
MOD = 1_000_000_007

class BerlekampMasseyKitamasa:
    def berlekampMassey(self, a):
        """返回最短递推系数 [c1, c2, ..., ck]"""
        pre_c = []; pre_i, pre_d = -1, 0; coef = []
        for i, v in enumerate(a):
            d = (v - sum(c * a[i-1-j] for j, c in enumerate(coef))) % MOD
            if d == 0: continue
            if pre_i < 0:
                coef = [0] * (i + 1); pre_i, pre_d = i, d; continue
            bias = i - pre_i
            old_len = len(coef)
            new_len = bias + len(pre_c)
            if new_len > old_len:
                tmp = coef[:]; coef += [0] * (new_len - old_len)
            delta = d * pow(pre_d, MOD - 2, MOD) % MOD
            coef[bias - 1] = (coef[bias - 1] + delta) % MOD
            for j, c in enumerate(pre_c):
                coef[bias + j] = (coef[bias + j] - delta * c) % MOD
            if new_len > old_len:
                pre_c = tmp; pre_i, pre_d = i, d
        return coef

    def kitamasa(self, coef, a, n):
        """根据递推系数和初始值，求第 n 项"""
        if n < len(a): return a[n] % MOD
        k = len(coef)
        def compose(ac, bc):
            res = [0] * k; temp = bc.copy()
            for c in ac:
                for j in range(k): res[j] = (res[j] + c * temp[j]) % MOD
                bk1 = temp[-1]
                for j in range(k-1, 0, -1): temp[j] = (temp[j-1] + bk1 * coef[j]) % MOD
                temp[0] = bk1 * coef[0] % MOD
            return res
        res_c = [0]*k; res_c[0] = 1
        c = [0]*k; c[1] = 1 if k > 1 else 0
        while n > 0:
            if n & 1: res_c = compose(c, res_c)
            c = compose(c, c.copy()); n >>= 1
        return sum(c * v for c, v in zip(res_c, a)) % MOD

bmk = BerlekampMasseyKitamasa()
fib = [0, 1, 1, 2, 3, 5, 8, 13]
coef = bmk.berlekampMassey(fib)
print("递推系数:", coef)       # [1, 1]
print("第 100 项:", bmk.kitamasa(coef, fib, 100))

总结 ​

算法时间复杂度适用场景
埃筛O(n log log n)快速筛质数
欧拉筛O(n)线性筛质数
质因数分解O(√n)分解质因数
快速幂O(log n)幂运算
矩阵快速幂O(k³ log n)线性递推
组合数预处理O(MAXN)质数模组合数
排列组合回溯O(n!) / O(C(n,k))枚举排列/组合
线性基O(n · 位数)最大异或和
幺半群滑动窗口O(1) 均摊滑窗聚合运算
根号分解O(n√n)按步长分治
几何类O(1) 操作点/直线/线段/圆
凸包O(n log n)求平面点集凸包
扫描线O(n² log n)矩形面积并
Berlekamp-MasseyO(k² log n)求最短递推式 + 第n项

模板源码：AlgorithmLearning