作业 3

若没有其他组员，答案为：none。

（若有组员，则列出每位成员的姓名与 SID。）

\(C\) 的次数为 \(\deg A + \deg B = 3 + 3 = 6\)，因此至少需要 \(6+1=7\) 个点才能通过逆 FFT 恢复 \(C\) 的系数。由于 \(n\) 必须是 2 的幂，故取 \(n = 8\)。

官方解答仅给出 \(n=8\) 而未列出 \(A\)、\(B\)、\(C\) 在各 8 次单位根处具体的取值。

我们证明 \(M_n(\omega^{-1})M_n(\omega)\) 的第 \((j,k)\) 个元素（行列均从 0 开始计数）在 \(j=k\) 时为 \(n\)，在 \(j\neq k\) 时为 \(0\)。

由矩阵乘法定义，第 \((j,k)\) 项为

\[\sum_{l=0}^{n-1} M_n(\omega^{-1})_{jl}\, M_n(\omega)_{lk} = \sum_{l=0}^{n-1} \omega^{-lj}\, \omega^{kl} = \sum_{l=0}^{n-1} \omega^{l(k-j)}.\]

若 \(j=k\)，则上式化为

\[\sum_{l=0}^{n-1} \omega^{0} = \sum_{l=0}^{n-1} 1 = n.\]

若 \(j\neq k\)，则 \(\omega^{k-j}\neq 1\)，由等比数列求和公式

\[\sum_{l=0}^{n-1} \omega^{l(k-j)} = \frac{1 - \omega^{n(k-j)}}{1 - \omega^{k-j}}.\]

由于 \(\omega\) 是 \(n\) 次单位根，对任意整数 \(m\) 都有 \(\omega^{nm}=1\)，故上式化为

\[\frac{1-1}{1-\omega^{k-j}} = 0.\]

另一种看法：因式分解 \(X^n - 1 = (X-1)(X^{n-1}+X^{n-2}+\dots+X+1)\)。\(\omega^{k-j}\) 是 \(X^n-1\) 的根；若 \(j\neq k\) 则 \(\omega^{k-j}\neq 1\)，故不是 \(X-1\) 的根，必为 \(X^{n-1}+\dots+1\) 的根，即

\[\omega^{(n-1)(k-j)} + \omega^{(n-2)(k-j)} + \dots + \omega^{k-j} + 1 = 0,\]

恰好就是要证的等式。

因此 \(\frac{1}{n}M_n(\omega^{-1})\,M_n(\omega) = I\)，即 \(\frac{1}{n}M_n(\omega^{-1})\) 是 \(M_n(\omega)\) 的逆矩阵。

观察得 \((M_n(\omega))^\dagger = M_n(\omega^{-1})\)（对每个元素取共轭即把 \(\omega\) 换成 \(\omega^{-1}=\bar{\omega}\)）。于是

\[\left(\frac{1}{\sqrt{n}}M_n(\omega)\right)\!\left(\frac{1}{\sqrt{n}}M_n(\omega)\right)^\dagger = \frac{1}{n}\,M_n(\omega)\,M_n(\omega^{-1}) = I,\]

其中最后一步用到 (b) 已证的 \(\frac{1}{n}M_n(\omega^{-1})M_n(\omega)=I\)（将 \(\omega\) 与 \(\omega^{-1}\) 角色对调即得）。故 \(\frac{1}{\sqrt{n}}M_n(\omega)\) 是酉矩阵。

设 \(c_0,\dots,c_{n-1}\) 为 \(C(x)\) 的系数。由课堂结论

\[M_n(\omega)\begin{bmatrix}c_0\\ c_1\\ \vdots \\ c_{n-1}\end{bmatrix} =\begin{bmatrix}C(1)\\ C(\omega)\\ \vdots \\ C(\omega^{n-1})\end{bmatrix}.\]

于是

\[\begin{bmatrix}c_0\\ c_1\\ \vdots \\ c_{n-1}\end{bmatrix} = M_n(\omega)^{-1}\begin{bmatrix}C(1)\\ C(\omega)\\ \vdots \\ C(\omega^{n-1})\end{bmatrix}.\]

由 (b) 知 \(M_n(\omega)^{-1} = \frac{1}{n}M_n(\omega^{-1})\)。因此，为求 \(C(x)\) 的系数，只需将 \(\frac{1}{n}M_n(\omega^{-1})\) 乘以取值向量 \([C(1)\; C(\omega)\; \dots\; C(\omega^{n-1})]^T\) 即可。

首先指出：若 \(\omega\) 是本原 \(n\) 次单位根，则 \(\omega^{-1}\) 也是本原 \(n\) 次单位根（因为 \(\omega^{-1}=e^{-i\theta}\) 当 \(\omega=e^{i\theta}\)）。我们在证明 \(M_n(\omega)\) 可拆成 4 块时只用到了「\(\omega\) 是本原 \(n\) 次单位根」这一性质，因此把 \(\omega\mapsto\omega^{-1}\) 代入仍成立，即得所要的重排后形式。

下面逐块直接验证（列按偶、奇重排）。

左上块：第 \((j,k)\) 个元素为 \(\omega^{-j\cdot(2k)} = \omega^{-2jk}\)。因 \(n\) 为偶数，\(\omega^2\) 是本原 \((n/2)\) 次单位根，故该元素正是 \(M_{n/2}(\omega^{-2})\) 的第 \((j,k)\) 项。

左下块：第 \((j,k)\) 个元素为 \(\omega^{-(n/2+j)(2k)} = \omega^{-nk}\,\omega^{-2jk}\)。由 \(\omega^n=1\) 得 \(\omega^{-nk}=1\)，故亦为 \(\omega^{-2jk}\)，与左上块相同。

右上块：第 \((j,k)\) 个元素为 \(\omega^{-j(2k+1)} = \omega^{-j}\,\omega^{-2jk}\)，即 \(\omega^{-j}\) 乘以左上块第 \((j,k)\) 项，故右上块 \(=\,\omega^{-j}M_{n/2}(\omega^{-2})\)。

右下块：第 \((j,k)\) 个元素为 \(\omega^{-(j+n/2)(2k+1)} = \omega^{-2jk}\,\omega^{-nk}\,\omega^{-j}\,\omega^{-n/2}\)。由 \(\omega^{-nk}=1\) 及 \((\omega^{-n/2})^2=1\) 但 \(\omega^{-n/2}\neq 1\) 可得 \(\omega^{-n/2}=-1\)，于是该元素 \(= -\omega^{-j}\,\omega^{-2jk}\)，即 \(-\omega^{-j}\) 乘以左上块对应项。

故重排后

\[M_n(\omega^{-1}) \rightarrow \begin{bmatrix} M_{n/2}(\omega^{-2}) & \omega^{-j}M_{n/2}(\omega^{-2}) \\ M_{n/2}(\omega^{-2}) & -\omega^{-j}M_{n/2}(\omega^{-2}) \end{bmatrix}.\]

附注：一般地，复数 \(x\) 与整数 \(n,m\) 下 \((x^n)^{1/m}\neq (x^{1/m})^n\)，故 \(x^{n/m}\) 并不良好定义。但此处 \(n/2\) 是整数，而整数次幂在复数上远比非整数次幂表现良好，故上述推导并无问题。

算法描述

构造两个向量：

• 向量 \(\vec{p} = [0, \dots, 0, b_0, \dots, b_{n-1}]\) 是「填充」向量，长度为 \(2n\)，由 \(n\) 个 0 后接 \(b\) 的每一位构成；

• 向量 \(\vec{f} = [(1-b_0), (1-b_1), \dots, (1-b_{n-1})]\) 是「翻转」向量，由 \(b\) 的每一位取反（0 变 1、1 变 0）构成。

对这两个向量运行跨相关算法。\(\vec{f}\) 与 \(\vec{p}[j : n+j]\) 的内积恰好是 \((n-j)\)-逆序的个数，对 \(j = 1, 2, \dots, n-1\) 报告即可。（注意第一个和最后一个内积会被忽略，因为不存在 0-逆序或 \(n\)-逆序。）

正确性证明

从下标 \(i\) 开始存在一个 \(k\)-逆序 当且仅当 \(b_i = 1\) 且 \(b_{i+k} = 0\)；这等价于 \(b_i(1-b_{i+k}) = 1\)。

考虑内积 \(\vec{f} \cdot \vec{p}[j : n+j] = \sum_{i=0}^{n-1} f_i p_{i+j}\)。由于当 \(i+j < n\) 时 \(p_{i+j} = 0\)，因此 \[\sum_{i=0}^{n-1} f_i p_{i+j} = \sum_{i=n-j}^{n-1} f_i p_{i+j} = \sum_{i=0}^{j-1} f_{i+n-j} p_{n+i}.\]

由 \(\vec{f}\) 与 \(\vec{p}\) 的定义，上式等于 \(\sum_{i=0}^{j-1} (1-b_{i+n-j}) b_i\)。根据第一句话，求和中的每一项在从 \(i\) 开始存在 \((n-j)\)-逆序时为 1，否则为 0，因此该和恰好是 \((n-j)\)-逆序的个数。

既然课上已知跨相关算法能正确计算这些内积，故整个算法正确。

运行时间

构造向量 \(\vec{p}\) 和 \(\vec{f}\) 共耗时 \(\Theta(n)\)；由课上结论，跨相关算法耗时 \(\Theta(n \log n)\)；整理输出结果耗时 \(\Theta(n)\)（或视实现为常数时间）。没有其他额外工作，故整体算法耗时 \(\Theta(n \log n)\)。

伪代码

CountKInversions(b, n):
    p ← array of length 2n, filled with 0
    for i ← 0 to n-1:
        p[n + i] ← b[i]
    f ← array of length n
    for i ← 0 to n-1:
        f[i] ← 1 - b[i]
    // 跨相关：等价于把 f 反转后与 p 做卷积
    corr ← CrossCorrelate(f, p)   // Θ(n log n)
    for j ← 1 to n-1:
        report corr[j] as number of (n-j)-inversions

复杂度：\(\Theta(n \log n)\) 时间，\(\Theta(n)\) 额外空间。

(a) \(O(nm)\) 算法

算法描述

对每个起始位置 \(i \in \{0, 1, \dots, m-n\}\)，检查子串 \(s[i : i+n-1]\) 是否与 \(g\) 匹配。每次检查耗时 \(O(n)\)，共有 \(O(m)\) 个起始位置，故总时间 \(O(mn)\)。

伪代码

NaiveMatch(g, n, s, m):
    matches ← empty list
    for i ← 0 to m - n:
        ok ← true
        for j ← 0 to n - 1:
            if g[j] ≠ s[i + j]:
                ok ← false; break
        if ok: append i to matches
    return matches

复杂度

时间 \(O(mn)\)，额外空间 \(O(1)\)（不计输出）。

---

(b) \(O(m \log m)\) 算法（比特串情形）

算法描述

令 \(g'\) 为把 \(g\) 中的 \(0\) 替换为 \(-1\) 得到的串；\(s'\) 为把 \(s\) 中的 \(0\) 替换为 \(-1\) 得到的串。

构造两个多项式：

• \(p_1(x) = a_0 + a_1 x + \cdots + a_{n-1} x^{n-1}\)，其中 \(a_d = g'(n-d-1)\)，即 \(g'\) 的系数顺序反转；

• \(p_2(x) = b_0 + b_1 x + \cdots + b_{m-1} x^{m-1}\)，其中 \(b_d = s'(d)\)。

令 \(p_3(x) = p_1(x) \cdot p_2(x) = c_0 + c_1 x + \cdots\)。则对任意 \(j \in \{0, 1, \dots, m-n\}\)，\(x^{n-1+j}\) 的系数为 \[c_{n-1+j} = \sum_{i=0}^{n-1} a_{n-1-i} b_{j+i} = \sum_{i=0}^{n-1} g'(i) s'(j+i),\] 恰好是 \(s'\) 中从下标 \(j\) 开始的子串与 \(g'\) 的内积。若两串在至多 \(k\) 个位置不同，则该内积至少为 \(n - 2k\)。

因此只需计算 \(p_3\)，并输出所有满足 \(c_{n-1+j} = n\) 的 \(j\)（此时子串与 \(g\) 完全匹配）。

正确性证明

在 \(\{-1, 1\}\) 编码下，\(g'(i) s'(j+i) = 1\) 当且仅当 \(g(i) = s(j+i)\)，否则为 \(-1\)。因此内积 \(\sum_{i=0}^{n-1} g'(i) s'(j+i)\) 等于（匹配位置数）减（不匹配位置数）。该和为 \(n\) 当且仅当所有 \(n\) 个位置都匹配，即子串与 \(g\) 完全相同。由 FFT 正确计算卷积系数，算法输出的恰好是所有匹配位置。

伪代码

BitstringMatch(g, n, s, m):
    g' ← array of length n, g'[i] ← (g[i] == 1 ? 1 : -1)
    s' ← array of length m, s'[i] ← (s[i] == 1 ? 1 : -1)
    // p1: 系数为 g' 的反转
    p1 ← array of length N (N 为不小于 m 的最小 2 的幂)
    for d ← 0 to n-1: p1[d] ← g'[n-1-d]
    // p2: 系数为 s'
    p2 ← array of length N
    for d ← 0 to m-1: p2[d] ← s'[d]
    q1 ← FFT(p1); q2 ← FFT(p2)           // O(m log m)
    for i ← 0 to N-1: q1[i] ← q1[i] * q2[i]   // 逐点乘, O(m)
    p3 ← InverseFFT(q1)                   // O(m log m)
    matches ← empty list
    for j ← 0 to m - n:
        if round(p3[n-1+j].real) == n:
            append j to matches
    return matches

复杂度

一次 FFT、一次逆 FFT 各 \(O(m \log m)\)，逐点乘与扫描系数各 \(O(m)\)，故总时间 \(O(m \log m)\)，空间 \(O(m)\)。

---

(c) \(O(m \log m)\) 算法（一般 \(\alpha\) 个字符情形）

算法描述

记 \(\omega_\alpha = e^{2\pi i / \alpha}\) 为 \(\alpha\) 次单位根。把字符编码为单位根的幂：

• \(g'(i) = \omega_\alpha^j\) 当 \(g(i) = c_j\)；

• \(s'(i) = \omega_\alpha^{-j}\) 当 \(s(i) = c_j\)。

于是 \(g'(i) s'(j+i) = 1\) 当且仅当 \(g(i) = s(j+i)\)，否则为一个不等于 1 的 \(\alpha\) 次单位根。

考虑内积 \(\sum_{i=0}^{n-1} g'(i) s'(j+i)\)：若子串与 \(g\) 完全匹配，则内积为 \(n\)；若存在任意不匹配位置，则和为 \(n\) 个 \(\alpha\) 次单位根之和，其中至少有一个不是 1。由于每个非 1 单位根的实部严格小于 1，\(n\) 个单位根之和为 \(n\) 当且仅当它们全是 1。因此内积为 \(n\) 当且仅当子串与 \(g\) 完全匹配。

按 (b) 的方式构造反转多项式 \(p_1\)（系数为 \(g'\) 反转）与 \(p_2\)（系数为 \(s'\)），用 FFT 计算 \(p_3 = p_1 \cdot p_2\)，输出所有满足 \(c_{n-1+j} = n\) 的 \(j\)。

伪代码

GeneralMatch(g, n, s, m, alpha, charToIdx):
    omega ← exp(2πi / alpha)
    g' ← array of length n
    for i ← 0 to n-1: g'[i] ← omega^(charToIdx[g[i]])
    s' ← array of length m
    for i ← 0 to m-1: s'[i] ← omega^(-charToIdx[s[i]])
    N ← 不小于 m 的最小 2 的幂
    p1 ← array of length N
    for d ← 0 to n-1: p1[d] ← g'[n-1-d]   // 反转
    p2 ← array of length N
    for d ← 0 to m-1: p2[d] ← s'[d]
    q1 ← FFT(p1); q2 ← FFT(p2)
    for i ← 0 to N-1: q1[i] ← q1[i] * q2[i]
    p3 ← InverseFFT(q1)
    matches ← empty list
    for j ← 0 to m - n:
        if |p3[n-1+j].real - n| < ε and |p3[n-1+j].imag| < ε:
            append j to matches
    return matches

复杂度

FFT 与逆 FFT 各 \(O(m \log m)\)，其余线性，故总时间 \(O(m \log m)\)，空间 \(O(m)\)。

主要思路：指数运算将乘法转化为加法。例如 \(x^3 \cdot x^2 = x^{2+3} = x^5\)。这提示我们将列表表示为多项式，其中列表的元素作为指数（这是一个非常常用的技巧，应当记住！）。例如，数组 \([1, 3]\) 可表示为 \(x^1 + x^3\)。将其与 \([2, 4]\) 对应的多项式 \(x^2 + x^4\) 相乘，得到 \(x^3 + 2x^5 + x^7\)。注意 \(3, 5, 7\) 分别对应 \([1,3]\) 中一个元素与 \([2,4]\) 中一个元素的所有可能和。

更正式地，定义

\[p(x) = x^{A[0]} + x^{A[1]} + \cdots + x^{A[n-1]}.\]

注意 \(p(x)^3\) 包含一系列项，每项形如 \(x^{A[i]} \cdot x^{A[j]} \cdot x^{A[k]} = x^{A[i]+A[j]+A[k]}\)。因此，我们只需检查 \(p(x)^3\) 是否包含 \(x^n\) 项。

正确性证明：观察

\[q(x) = p(x)^3 = \left(\sum_{0 \leq i < n} x^{A[i]}\right)^3 = \sum_{0 \leq i,j,k < n} x^{A[i]+A[j]+A[k]}.\]

因此，\(q\) 中 \(x^n\) 的系数非零当且仅当存在下标 \(i, j, k\) 使得 \(A[i] + A[j] + A[k] = n\)。故算法正确。（实际上，\(x^n\) 的系数告诉我们满足条件的三元组 \((i,j,k)\) 有多少个。）

复杂度分析：构造 \(p(x)\) 显然耗时 \(O(n)\)。\(p(x)\) 是次数至多 \(n = O(n)\) 的多项式。因此用 FFT 做两次乘法计算 \(q(x)\) 耗时 \(O(n \log n)\)。最后查找 \(x^n\) 的系数耗时 \(O(1)\)，故总复杂度为 \(O(n \log n)\)。

备注：本题保证数组元素在 \(0 \ldots n\) 范围内。如果没有这个保证，上述 FFT 方法会变得低效（因为多项式次数与 \(A\) 中最大元素同阶）。此时容易得到 \(O(n^2)\) 算法，但尚未发现更快的算法。特别是，是否存在比 \(O(n^2)\) 更快的算法是一个著名的开放问题（称为 3SUM 问题）。该问题被广泛研究，因为它与计算几何中的许多问题密切相关。

Q1.1) 单位根

def roots_of_unities(n):
    principal_root = np.exp(2 * np.pi * 1j / n)
    roots = principal_root ** np.arange(n)
    return roots

替代解法（求 \(x^n - 1 = 0\) 的根并按辐角排序）：

def roots_of_unities(n):
    poly = np.zeros(n + 1)
    poly[0] = 1
    poly[-1] = -1
    roots = np.roots(poly)
    sorted_roots = sorted(roots, key=lambda x: np.angle(x) % (2 * np.pi))
    return sorted_roots

Q1.2) 朴素 DFT

def dft_naive(coeffs):
    n = hyperceil(len(coeffs))
    p = np.poly1d(coeffs[::-1])
    return [p(omega) for omega in roots_of_unities(n)]

Q1.3) FFT

def fft(coeffs, roots):
    n = len(coeffs)
    assert n == hyperceil(n)
    if n == 1:
        return coeffs[:]
    e_coeff = coeffs[0::2]
    o_coeff = coeffs[1::2]
    e = fft(e_coeff, roots[::2])
    o = fft(o_coeff, roots[::2])
    res = np.zeros(n, dtype=np.complex_)
    for i in range(n // 2):
        res[i] = e[i] + roots[i] * o[i]
        res[i + n // 2] = e[i] + roots[i + n // 2] * o[i]
    return res

Q1.4) 逆 FFT

def ifft(vals, roots):
    n = len(vals)
    assert n == hyperceil(n)
    iroots = np.array(roots) ** -1
    return np.array(fft(vals, iroots) / n)

Q1.5) 多项式乘法 (可选)

def pad(coeffs, to):
    return np.pad(coeffs, (0, to - len(coeffs)), "constant", constant_values=0)

def poly_multiply(coeffs1, coeffs2):
    n = hyperceil(len(coeffs1) + len(coeffs2) - 1)
    v1 = fft(pad(coeffs1, n), roots_of_unities(n))
    v2 = fft(pad(coeffs2, n), roots_of_unities(n))
    v = v1 * v2
    return ifft(v, roots_of_unities(n))

注意：FFT 使用复数运算会导致浮点误差，整数多项式乘法结果需取整并去除尾部零：

def round_complex_to_int(lst):
    return [round(x.real) for x in lst]

def zero_pop(lst):
    return np.trim_zeros(lst, "b")

Q2) 预测评估

def evaluate_predictions(coeffs, labels):
    coeffs_shift = np.insert(coeffs, 0, 0)
    n = hyperceil(len(coeffs_shift) + len(labels) - 1)
    C = fft(pad(coeffs_shift, n), roots_of_unities(n))
    Y = fft(pad(labels, n), roots_of_unities(n))
    P = C * Y
    p = ifft(P, roots_of_unities(n))[1 : len(labels)]
    return np.mean((p - labels[1:]) ** 2)

算法思路：预测 \(p_t = \sum_{k=0}^{t-1} c_k y_{t-1-k}\) 是系数序列 \(c\)（前插 0 移位后）与标签序列 \(y\) 的卷积。通过 FFT 将两序列转为点值形式，逐点相乘后再逆 FFT 得到所有 \(p_t\)，最后线性扫描计算均方误差。总复杂度 \(O(n \log n)\)。