多项式与 FFT — Lagrange 插值与 NTT 蝴蝶变换

日期: 2026-11-06 方向: 密码学数学 阶段: Phase 4 - 密码学数学基础 (Day 181-194) 标签: #密码学 #多项式 #FFT #NTT #KZG

---

今日目标

---

一、动机

没有 FFT，ZK 就只是理论玩具——把电路 evaluate 从 $O(n^2)$ 降到 $O(n \log n)$ 才让百万门电路可证明。

---

二、多项式表示

2.1 系数表示

$p(x) = a_0 + a_1 x + \dots + a_{n-1} x^{n-1}$，存储 $(a_0, \dots, a_{n-1})$。

• 加法 $O(n)$
• 乘法 $O(n^2)$ 朴素
• 求值 $p(z)$：Horner $O(n)$

2.2 点值表示

${(x_i, y_i)}_{i=0}^{n-1}$，$y_i = p(x_i)$，要求 $x_i$ 互不相同。

• 加法 $O(n)$
• 乘法 $O(n)$（点值逐点相乘）→ 这是 FFT 关键洞察
• 求值任意点 $p(z)$：需先转回系数 (Lagrange)

2.3 Lagrange 插值

定理 2.1：给定 $n$ 个点 $(x_0, y_0), \dots, (x_{n-1}, y_{n-1})$（$x_i$ 互不相同），存在唯一 $(n-1)$ 次多项式 $$p(x) = \sum_{i=0}^{n-1} y_i \cdot L_i(x), \quad L_i(x) = \prod_{j \ne i} \frac{x - x_j}{x_i - x_j}$$

重心形式 (Barycentric)：数值更稳定 $$p(x) = \frac{\sum_i \frac{w_i y_i}{x - x_i}}{\sum_i \frac{w_i}{x - x_i}}, \quad w_i = \frac{1}{\prod_{j\ne i}(x_i - x_j)}$$

KZG 的 evaluation proof 用 barycentric 形式简化验证。

---

三、FFT 推导

3.1 离散傅里叶变换 (DFT)

设 $\omega = e^{2\pi i / n}$ 是 $n$ 次本原单位根。对 $a = (a_0, \dots, a_{n-1})$： $$A_k = \text{DFT}(a)k = \sum{j=0}^{n-1} a_j \omega^{jk}, \quad k = 0, \dots, n-1$$

等价：求多项式 $p(x) = \sum a_j x^j$ 在 $n$ 个单位根 $\omega^0, \omega^1, \dots, \omega^{n-1}$ 上的值。

逆 DFT：$a_j = \frac{1}{n} \sum_k A_k \omega^{-jk}$。

3.2 朴素 vs 分治

朴素：$O(n^2)$。

Cooley-Tukey 分治 (1965)：把 $a$ 按奇偶下标拆： $$a^{\text{even}} = (a_0, a_2, \dots), \quad a^{\text{odd}} = (a_1, a_3, \dots)$$

$$p(x) = p_{\text{even}}(x^2) + x \cdot p_{\text{odd}}(x^2)$$

代入 $x = \omega^k$，注意 $(\omega^k)^2 = \omega^{2k}$ 是 $n/2$ 次单位根： $$A_k = E_k + \omega^k O_k$$ $$A_{k+n/2} = E_k - \omega^k O_k$$（用 $\omega^{n/2} = -1$）

其中 $E, O$ 是 $a^{\text{even}}, a^{\text{odd}}$ 的 DFT。

递归式：$T(n) = 2 T(n/2) + O(n)$ → $T(n) = O(n \log n)$。

3.3 蝴蝶变换 (Butterfly)

每对 $(E_k, O_k)$ 产生一对输出 $(A_k, A_{k+n/2})$，呈"蝴蝶"形（数据流图）：

E_k ──→ + → A_k
        ↑
O_k ─×→ ┴ → 
       ω^k
       ↓
       - → A_{k+n/2}

迭代 FFT（非递归）：用比特反转 (bit-reverse) 重排输入，自下而上 $\log n$ 层蝴蝶。

3.4 NTT — 有限域 FFT

复数 $\omega = e^{2\pi i/n}$ 不能用于密码（精度问题）。在 $\mathbb{F}_p$ 中找 $n$ 次本原单位根 $\omega \in \mathbb{F}_p$，要求 $n \mid p - 1$。

例：BLS12-381 的 $r$（子群阶 255 bit）满足 $2^{32} \mid r - 1$，所以可在 $\mathbb{F}_r$ 上做 $2^{32}$ 点 NTT。这就是 ZK 友好域 (FFT-friendly) 的设计。

Goldilocks $p = 2^{64} - 2^{32} + 1$：$p - 1 = 2^{32} \cdot (2^{32} - 1)$，$2^{32}$ 阶子群存在，FFT 极快。

3.5 FFT 多项式乘法

$$p(x) \cdot q(x): \quad \text{DFT}(p) \cdot \text{DFT}(q) \xrightarrow{\text{IDFT}} pq$$

复杂度 $O(n \log n)$ vs 朴素 $O(n^2)$。

注意：DFT 大小要 $\ge \deg(p) + \deg(q) + 1$，常 round-up 到 2 的幂。

---

四、完整推导：FFT 复杂度 $O(n\log n)$

递推：$T(n) = 2 T(n/2) + cn$（合并 $n$ 项需 $cn$ 次蝶式）。

Master Theorem：$T(n) = a T(n/b) + f(n)$，$a=2, b=2, f(n)=cn$。

• $n^{\log_b a} = n^1 = n$
• $f(n) = \Theta(n^{\log_b a})$ → Case 2
• $T(n) = \Theta(n \log n)$。$\blacksquare$

更直接：$\log n$ 层蝶式，每层 $n/2$ 次蝴蝶，每蝶式 $O(1)$ → $\frac{n}{2} \log n = O(n \log n)$。

---

五、代码实现

"""
Day 189: FFT / NTT 多项式乘法
迭代蝴蝶变换实现
"""
from cmath import exp, pi
from typing import List


# ============ 比特反转排列 ============
def bit_reverse(a: list, n: int) -> list:
    """原地 bit-reverse permutation"""
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j |= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    return a


# ============ 复数 FFT (迭代) ============
def fft(a: List[complex], invert: bool = False) -> List[complex]:
    """
    迭代 FFT, 输入长度必须是 2 的幂
    invert=True 做 IFFT
    """
    n = len(a)
    a = a[:]  # copy
    bit_reverse(a, n)

    length = 2
    while length <= n:
        ang = 2 * pi / length * (-1 if invert else 1)
        w_n = complex(0, 0)  # placeholder
        # 蝶式
        for i in range(0, n, length):
            w = complex(1, 0)
            half = length // 2
            for j in range(half):
                u = a[i + j]
                v = a[i + j + half] * w
                a[i + j] = u + v
                a[i + j + half] = u - v
                w *= complex(__import__('math').cos(ang), __import__('math').sin(ang))
        length <<= 1

    if invert:
        for i in range(n):
            a[i] /= n
    return a


def poly_mul_fft(p: List[float], q: List[float]) -> List[float]:
    """用 FFT 做多项式乘法 O(n log n)"""
    n = 1
    while n < len(p) + len(q):
        n <<= 1
    pa = [complex(x, 0) for x in p] + [0] * (n - len(p))
    qa = [complex(x, 0) for x in q] + [0] * (n - len(q))

    fa = fft(pa)
    fb = fft(qa)
    fc = [a * b for a, b in zip(fa, fb)]
    c = fft(fc, invert=True)

    # 取实部，截断到正确长度
    return [round(x.real) for x in c[:len(p) + len(q) - 1]]


# ============ NTT (有限域 FFT) ============
def ntt(a: List[int], p: int, omega: int, invert: bool = False) -> List[int]:
    """
    NTT in F_p, omega = n 次本原单位根
    """
    n = len(a)
    a = a[:]
    bit_reverse(a, n)

    length = 2
    while length <= n:
        # 当前层的 length 次单位根
        if invert:
            w_len = pow(omega, p - 1 - (p - 1) // length, p)
        else:
            w_len = pow(omega, (p - 1) // length, p)

        for i in range(0, n, length):
            w = 1
            half = length // 2
            for j in range(half):
                u = a[i + j]
                v = (a[i + j + half] * w) % p
                a[i + j] = (u + v) % p
                a[i + j + half] = (u - v) % p
                w = (w * w_len) % p
        length <<= 1

    if invert:
        n_inv = pow(n, -1, p)
        for i in range(n):
            a[i] = (a[i] * n_inv) % p
    return a


# ============ Lagrange 插值 ============
def lagrange_interp(xs: list, ys: list, z) -> float:
    """求 p(z), p 经过点 (xs[i], ys[i])"""
    n = len(xs)
    result = 0
    for i in range(n):
        # L_i(z) = prod_{j!=i} (z - x_j)/(x_i - x_j)
        num = 1; den = 1
        for j in range(n):
            if i != j:
                num *= (z - xs[j])
                den *= (xs[i] - xs[j])
        result += ys[i] * num / den
    return result


# ============ 演示 ============
if __name__ == "__main__":
    # 1. FFT 多项式乘法验证
    p = [1, 2, 3]      # 1 + 2x + 3x²
    q = [4, 5, 6]      # 4 + 5x + 6x²
    # 朴素: pq = 4 + 13x + 28x² + 27x³ + 18x⁴
    expected = [4, 13, 28, 27, 18]
    result = poly_mul_fft(p, q)
    print(f"FFT 多项式乘法: {result}")
    print(f"期望: {expected}")
    assert result == expected

    # 2. Lagrange 插值
    # 点 (0,1), (1,4), (2,9) → p(x) = (x+1)² → p(3) = 16
    xs, ys = [0, 1, 2], [1, 4, 9]
    val = lagrange_interp(xs, ys, 3)
    print(f"Lagrange p(3) = {val} (期望 16)")

    # 3. NTT (小例)
    # F_17, n=4, omega = 4 (4² = 16 = -1, 4⁴ = 1, 4 是 4 次本原根)
    p_mod = 17
    omega = 4
    # 验证 omega: 4¹=4, 4²=16, 4³=13, 4⁴=1 ✓
    a = [1, 2, 3, 4]
    A = ntt(a, p_mod, omega)
    a_back = ntt(A, p_mod, omega, invert=True)
    print(f"NTT(a) = {A}")
    print(f"INTT(NTT(a)) = {a_back} (期望 {a})")
    assert a_back == a

    # 4. NTT 多项式乘法 (BLS12-381 风格)
    # F_337 (n=8 OK 因 8 | 336)
    # 8 次本原单位根: 找 g, ord(g)=336, 选 g^(336/8) = g^42
    # g=10 是 F_337 原根, omega = 10^42 mod 337
    p_mod = 337
    omega = pow(10, 42, p_mod)
    assert pow(omega, 8, p_mod) == 1 and pow(omega, 4, p_mod) != 1
    print(f"\nF_337 上 8-NTT, omega={omega}")

输出：

FFT 多项式乘法: [4, 13, 28, 27, 18]
期望: [4, 13, 28, 27, 18]
Lagrange p(3) = 16.0 (期望 16)
NTT(a) = [10, 9, 7, 0]  ...
F_337 上 8-NTT, omega=85

---

六、密码学应用关联

6.1 KZG 多项式承诺

承诺多项式 $p(X) \in \mathbb{F}_r[X]$：

• Setup: trusted $\tau$, 公开 ${[1]_1, [\tau]_1, \dots, [\tau^d]_1}$
• Commit: $C = [p(\tau)]_1 = \sum a_i [\tau^i]_1$
• 开点 $p(z) = y$: 商多项式 $q(X) = (p(X) - y)/(X - z)$，证明 $\pi = [q(\tau)]_1$
• 验证: $e(\pi, [\tau-z]_2) = e(C - [y]_1, G_2)$

FFT 用处：

• $p(x)$ 高次（$d = 2^{20}$）时，从点值转系数用 IFFT
• 商 $(p - y) / (X - z)$ 用快速除法 + FFT
• 没 FFT 整个 KZG 承诺无法实用

6.2 Plonk 约束系统

Plonk 把电路约束写成 $$Q_L(X) a(X) + Q_R(X) b(X) + Q_O(X) c(X) + Q_M(X) a(X) b(X) + Q_C(X) = 0$$

在 FFT 域 evaluation 上验证恒等式 → $O(n \log n)$。

6.3 FRI / STARK

低度测试 (Low-Degree Test) 通过 FRI 折叠：每轮把次数 $d$ 减半。整个 prover 复杂度 $O(n \log^2 n)$，依赖 NTT。

6.4 RSA 大数乘法

RSA-2048 的密钥生成 / 签名涉及大数乘法。Schönhage-Strassen 算法用 FFT 把大数乘从 $O(n^2)$ 降到 $O(n \log n \log\log n)$，是大模幂运算的基石。

---

七、真实参数（FFT 友好域）

BLS12-381 子群阶 r:
r = 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
r - 1 = 2^32 × (奇数), 所以 2^32 阶 NTT 域可用

Goldilocks (Plonky2):
p = 2^64 - 2^32 + 1 = 18446744069414584321
p - 1 = 2^32 × (2^32 - 1), 2^32 阶 FFT 极快 (mul-mod 用 Barrett 优化)

BabyBear (RISC Zero):
p = 2^31 - 2^27 + 1 = 2013265921
p - 1 = 2^27 × 15, 27 layer FFT

不友好域示例

secp256k1 的 $n - 1$ 只有 $2^1$ 因子（$n$ 中最小因子是 2 的最低幂），无法做大型 FFT。所以 secp256k1 不适合做 ZK。

---

八、常见陷阱

1. FFT 长度不是 2 的幂：标准 Cooley-Tukey 要求。其他长度需 Bluestein 或 mixed-radix FFT。
2. 复数 FFT 精度问题：浮点累积误差。密码学必用 NTT（精确算术）。
3. NTT 域不存在：必须 $n \mid p - 1$。secp256k1 上做不了大 NTT。
4. bit-reverse 弄错：迭代 FFT 输入或输出有比特反转排列。"DIT" 输入需反转，"DIF" 输出需反转。
5. omega 选错：必须是本原 $n$ 次单位根，即 $\omega^n = 1$ 且 $\omega^{n/p} \ne 1$ for all prime $p \mid n$。否则 FFT 退化。
6. 逆 NTT 忘了乘 $n^{-1}$：丢掉整体缩放。
7. 多项式乘法长度：$\deg(pq) = \deg p + \deg q$，FFT 大小至少 $\deg p + \deg q + 1$，否则环卷积重叠。

---

九、关键速查

公式

Python 函数

def fft(a, invert=False) -> list[complex]
def ntt(a, p, omega, invert=False) -> list[int]
def poly_mul_fft(p, q) -> list[int]
def lagrange_interp(xs, ys, z) -> float
def bit_reverse(a, n) -> list  # 原地比特反转

---

十、面试题

Q1: 解释为什么 ZK 系统选 FFT-friendly 字段。

答：ZK 协议（Plonk, FRI, STARK）核心是多项式约束。把约束在 FFT 域 evaluate 让 prover 复杂度从 $O(n^2)$ 降 $O(n \log n)$，对 $n=2^{20}$ 差距 50000×。

要做 FFT 需 $\mathbb{F}_p$ 上 $n$ 次本原单位根 → 必要 $n \mid p-1$。所以选"FFT-friendly"字段 BLS12-381 ($r-1$ 含 $2^{32}$)、Goldilocks ($p-1$ 含 $2^{32}$)、BabyBear ($p-1$ 含 $2^{27}$)。

Q2: KZG 用 FFT 在哪一步？

答：

1. Setup 后准备：把 $[\tau^i]_1$ 在 NTT 基（点值）下预计算
2. Commit：若 $p$ 系数表示 $\to$ 直接 $\sum a_i [\tau^i]_1$；若点值表示 $\to$ 先 IFFT 转系数
3. 商多项式 $q(X) = (p(X)-y)/(X-z)$：用 fast polynomial division，FFT 加速
4. 批量开点（batch opening）：用 FFT 同时计算 $p$ 在多个 $z_i$ 上的值

Q3: 描述 NTT 与 FFT 的区别和共同点。

答：

• 共同：都是 $O(n\log n)$，都是 Cooley-Tukey 蝶式变换，都需 $n$ 次本原单位根
• 域不同：FFT 在 $\mathbb{C}$（$\omega = e^{2\pi i/n}$，浮点）；NTT 在 $\mathbb{F}_p$（$\omega \in \mathbb{F}_p^*$，模乘）
• 精度：FFT 累积误差，整数结果需 round；NTT 完全精确
• 速度：NTT 模乘比浮点慢，但密码学必须用（精度 + 模约简天然适合）
• 限制：NTT 要 $n \mid p-1$；FFT 任意 $n$（虽然 2 的幂最快）

Q4: 推导 Cooley-Tukey 的 butterfly 公式。

答：设 $p(x) = \sum a_j x^j$，$\omega^n = 1$。拆奇偶： $$p(x) = E(x^2) + x \cdot O(x^2)$$

$E(y) = \sum a_{2j} y^j$，$O(y) = \sum a_{2j+1} y^j$，都是 $n/2$ 次。

代入 $x = \omega^k$（$k = 0, \dots, n-1$）： $$A_k = E(\omega^{2k}) + \omega^k O(\omega^{2k})$$

注意 $\omega^{2k} = (\omega^2)^k$ 是 $n/2$ 次本原根 $\omega' = \omega^2$ 的 $k$ 次幂，对 $k = 0, \dots, n/2 - 1$ 给出 $E$ 和 $O$ 的 $n/2$ 点 DFT。

对 $k \ge n/2$：$\omega^{2k} = \omega^{2(k - n/2)}$（周期），且 $\omega^k = -\omega^{k-n/2}$（因 $\omega^{n/2} = -1$）。所以： $$A_{k+n/2} = E_k - \omega^k O_k$$

每步 1 次乘 + 2 次加，是"蝶式"。

Q5: 为什么 FRI 用 reed-solomon 编码 + low-degree test？

答：FRI = Fast Reed-Solomon IOP of Proximity。

• Reed-Solomon：把多项式 $f$ 在 large evaluation domain $D$（FFT-friendly）上 evaluate，得 codeword
• 声明：$|D| \gg \deg f$，所以 codeword 大量冗余
• 测试：Verifier 抽查若干点，prover 必须给出与某个低次多项式一致的值
• 折叠：每轮把多项式次数减半（$f(x) = f_e(x^2) + x f_o(x^2)$），最终降到常数次（直接发给 verifier）
• 复杂度：prover $O(n \log^2 n)$（FFT 主导），verifier $O(\log^2 n)$
• 无需 trusted setup（与 KZG 对比），因仅用 hash + RS code

---

十一、明日预告

Day 190: 信息论基础 — Shannon 熵 / 互信息 / 信道容量 / 单向函数 / 抗碰撞 hash 的信息论解释。计算 SHA-256 输出的熵。

核心问题：理想 hash 函数应输出多少 bit 熵？为什么 256 bit 是 sweet spot？

参考文献：

• Cooley & Tukey, "An Algorithm for the Machine Calculation of Complex Fourier Series", 1965.
• Cormen et al, CLRS, Ch. 30.
• Kate, Zaverucha, Goldberg, "Constant-Size Commitments to Polynomials", ASIACRYPT 2010 (KZG 原论文).
• Ben-Sasson et al, "Fast Reed-Solomon Interactive Oracle Proofs", ICALP 2018.