斐波那契：矩阵快速幂、通项与扩域

平凡的，我们有斐波那契数列：

F_0 = 0,\quad F_1 = 1,\quad F_{n+1} = F_n + F_{n-1}\ (n \ge 1)

递推的矩阵形式

\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix}

记

A = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}

则

\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = A^n \begin{pmatrix} F_1 \\ F_0 \end{pmatrix} = A^n \begin{pmatrix} 1 \\ 0 \end{pmatrix}

因此算 $F_n$ 等价于算 $A^n$ 。

复杂度模型

我们约定：

设 $M(t)$ 为两个 $t$ 位整数相乘的代价； $M$ 单调非降。
$F_n \sim \lambda_1^n / \sqrt{5}$ （ $\lambda_1$ 是其中一个特征值，证明过程见后文），故 $F_n$ 的二进制位数

N := \bigl\lfloor\log_2 F_n\bigr\rfloor + 1 = n\log_2\lambda_1 + O(1) = \Theta(n).

二进制快速幂算 $g^n$ （按最优左→右计数）：恰做 $\lfloor\log_2 n\rfloor$ 次平方，以及 $\mathrm{popcount}(n)-1$ 次一般乘，环乘法次数为

\mu(n) := \lfloor\log_2 n\rfloor + \mathrm{popcount}(n) - 1.

于是

\lfloor\log_2 n\rfloor \le \mu(n) \le 2\lfloor\log_2 n\rfloor,

即 $\mu(n) = \Theta(\log n)$ 。

矩阵快速幂时，中间矩阵的各元素的位数都不超过 $F_{n+1}$ ，与 $F_n$ 同阶位数，故其中的大整数乘按 $M(N)$ 计代价。
扩域法的中间量是 $Y_n = F_n\cdot 2^{n-1}$ ，位数

N' := N + (n-1) = n(1+\log_2\lambda_1) + O(1) = \Theta(n),

大整数乘按 $M(N')$ 计——与 $M(N)$ 同阶，但常数更大，位数比 $N'/N\to 2.44$ （详见 §位数比的解析效应）。

1. 矩阵快速幂

对 $A$ 做二进制快速幂得 $A^n$ ，再右乘 $(1,0)^\top$ ，读出第二分量即 $F_n$ 。

$2\times 2$ 矩阵朴素乘法恰含 $8$ 次 元素乘法（ $4$ 个点积，各 $2$ 次乘）。故整数乘次数为 $8\,\mu(n)$ ，总代价

\Theta\bigl(8\,M(N)\,\mu(n)\bigr) = \Theta\bigl(M(N)\log n\bigr).

这是最直接的精确 $O(\log n)$ 次环运算算法。

2. 特征值与对角化

目标是把 $A$ 对角化成 $A = PDP^{-1}$ ，得到

A^n = (PDP^{-1})^n = PD^n P^{-1}.

因为对角阵的幂次有一条极优雅的性质： $D = \mathrm{diag}(a,b)$ 时 $D^n = \mathrm{diag}(a^n, b^n)$ 。这可以极大地减轻运算负担。

特征方程

A\vec{x} = \lambda\vec{x} \quad\Rightarrow\quad (A - \lambda I)\vec{x} = 0 \quad\Rightarrow\quad \det(A - \lambda I) = 0

\begin{vmatrix} 1-\lambda & 1 \\ 1 & -\lambda \end{vmatrix} = (1-\lambda)(-\lambda) - 1 = 0

即

\lambda^2 - \lambda - 1 = 0.

解得两个特征值

\lambda_1 = \frac{1+\sqrt{5}}{2},\qquad \lambda_2 = \frac{1-\sqrt{5}}{2}.

易见 $\lambda_1 > 1$ ， $|\lambda_2| < 1$ ，且 $\lambda_1 + \lambda_2 = 1$ ， $\lambda_1\lambda_2 = -1$ ， $\lambda_1 - \lambda_2 = \sqrt{5}$ 。

特征向量

由 $(A - \lambda I)\vec{x} = 0$ ，第二行给出 $x_1 = \lambda x_2$ ，取特解

\vec{x}_1 = \begin{pmatrix} \lambda_1 \\ 1 \end{pmatrix},\quad \vec{x}_2 = \begin{pmatrix} \lambda_2 \\ 1 \end{pmatrix}.

令 $P$ 的列为这两个特征向量：

P = \begin{pmatrix} \lambda_1 & \lambda_2 \\ 1 & 1 \end{pmatrix},\quad D = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix}.

则 $AP = PD$ ，即 $A = PDP^{-1}$ 。

求逆并推出通项

\det(P) = \lambda_1 - \lambda_2 = \sqrt{5},

P^{-1} = \frac{1}{\sqrt{5}} \begin{pmatrix} 1 & -\lambda_2 \\ -1 & \lambda_1 \end{pmatrix}.

于是

\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = P D^n P^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix}.

先算

P^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{5}} \begin{pmatrix} 1 \\ -1 \end{pmatrix},\quad D^n P^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{5}} \begin{pmatrix} \lambda_1^n \\ -\lambda_2^n \end{pmatrix},

再左乘 $P$ ：

\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \frac{1}{\sqrt{5}} \begin{pmatrix} \lambda_1 & \lambda_2 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} \lambda_1^n \\ -\lambda_2^n \end{pmatrix} = \frac{1}{\sqrt{5}} \begin{pmatrix} \lambda_1^{n+1} - \lambda_2^{n+1} \\ \lambda_1^n - \lambda_2^n \end{pmatrix}.

得到 Binet 公式（对一切整数 $n\ge 0$ 精确成立）：

F_n = \frac{\lambda_1^n - \lambda_2^n}{\sqrt{5}} = \frac{1}{\sqrt{5}} \left( \left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right).

浮点快速幂

在浮点语义下，对 $\lambda_1^n$ 做快速幂再除以 $\sqrt{5}$ 、四舍五入，环运算次数仍是 $\mu(n)=\Theta(\log n)$ ，每次是机器浮点乘——这是有明确适用范围的精确算法：当浮点算出的 $\widetilde{x}$ 仍落在 $(F_n-1/2,\ F_n+1/2)$ 内时，四舍五入就得到 $F_n$ 。

能力边界由精度决定。IEEE-754 binary64 有效位数约 $53$ bit； $\lambda_1^n / \sqrt{5}$ 的量级约为 $F_n$ ，相对误差一旦大到使 $\widetilde{x}$ 跨过半整数边界，四舍五入就会错。后文实测：对本机 f64 快速幂实现，与精确值一致到 $n \le 75$ ，从 $n=76$ 起失准。

因此：小 $n$ 用浮点 Binet 最省事；要任意大的精确 $F_n$ ，必须离开浮点，改做整数算术。下一节正是通项在整数环上的落地。

3. 扩域

已知 $\forall n,\ F_n \in \mathbb{Z}$ ，Binet 里的 $\sqrt{5}$ 最终消掉。把计算放进环 $\mathbb{Z}[\sqrt{5}]$ ，全程整数，就不受浮点精度限制。

设

d_1 = 1 + \sqrt{5},\quad d_2 = 1 - \sqrt{5},

则 $\lambda_1 = d_1/2$ ， $\lambda_2 = d_2/2$ ，且

F_n = \frac{1}{\sqrt{5}} \left( \left(\frac{d_1}{2}\right)^n - \left(\frac{d_2}{2}\right)^n \right).

运算法则类似复数：把 $\sqrt{5}$ 当作「虚部」单位（ $\sqrt{5}^2 = 5$ ）。

$d_1$ 、 $d_2$ 共轭，故类似地，存在 $X_n, Y_n \in \mathbb{Z}$ 使

d_1^n = X_n + Y_n\sqrt{5},\qquad d_2^n = X_n - Y_n\sqrt{5}.

于是

\lambda_1^n - \lambda_2^n = \frac{d_1^n - d_2^n}{2^n} = \frac{2 Y_n \sqrt{5}}{2^n} = \frac{Y_n \sqrt{5}}{2^{n-1}},

F_n = \frac{\lambda_1^n - \lambda_2^n}{\sqrt{5}} = \frac{Y_n}{2^{n-1}}.

算法：快速幂算 $d_1^n = X_n + Y_n\sqrt{5}$ ，取出系数 $Y_n$ ，再除以 $2^{n-1}$ 。因 $F_n \in \mathbb{Z}$ ， $Y_n$ 必被 $2^{n-1}$ 整除；实现时用整数除法。

拓展：环乘法的整数乘次数

设 $\alpha = a + b\sqrt{5}$ ， $\beta = c + d\sqrt{5}$ 。朴素乘法

\alpha\beta = (ac + 5bd) + (ad + bc)\sqrt{5}

恰需 $4$ 次整数乘。但是，我们发现用一次分配律可以利用上之前计算的 $ac$ 和 $bd$ ，这样可以减少一次乘法，代价是两次加法：

ad + bc = (a+b)(c+d) - ac - bd,

得到

\alpha\beta = \bigl(ac + 5bd\bigr) + \bigl((a+b)(c+d) - ac - bd\bigr)\sqrt{5},

恰需 $3$ 次大整数乘（ $ac$ 、 $bd$ 、 $(a+b)(c+d)$ ）；其中 $\times 5$ 用移位加法 $5v=(v\ll 2)+v$ ，不计入大整数乘。

平方更简单：

(a + b\sqrt{5})^2 = (a^2 + 5b^2) + 2ab\sqrt{5},

也是 $3$ 次（ $a^2$ 、 $b^2$ 、 $ab$ ）。

中间量位数是 $N'=\log_2(Y_n)=\log_2(F_n\cdot 2^{n-1})=N+n-1$ ，故总代价

\Theta\bigl(3\,M(N')\,\mu(n)\bigr) = \Theta\bigl(M(N')\log n\bigr).

与矩阵的 $\Theta\bigl(8\,M(N)\,\mu(n)\bigr)$ 同属 $\Theta(M(n)\log n)$ 。单看乘法次数时 $3<8$ ，但 $N'>N$ ，乘次数优势会被更大的中间位数吃掉——这跟乘法的效率息息相关，具体分析见下。

位数比的解析效应

两种算法的总代价分别是 $8\,\mu(n)\,M(N)$ 与 $3\,\mu(n)\,M(N')$ 。环乘次数 $\mu(n)$ 在比值中消去，关键是 $M(N')/M(N)$ 。假设乘法代价是 $\alpha$ -次多项式增长：

M(t) = \Theta(t^\alpha),\qquad \alpha \ge 1.

则位数比的效应可以解析地写出来：

\frac{M(N')}{M(N)} = \left(\frac{N'}{N}\right)^\alpha,\qquad \frac{T_{\mathrm{ext}}}{T_{\mathrm{matrix}}} = \frac{3}{8}\left(\frac{N'}{N}\right)^\alpha.

$n\to\infty$ 时求位数比的极限（ $N=\log_2(F_n)+O(1)=\log_2(\lambda_1^n)+O(1)=n\log_2\lambda_1+O(1)$ ， $N'=n(1+\log_2\lambda_1)+O(1)$ ）：

\frac{N'}{N} = \frac{n(1+\log_2\lambda_1)+O(1)}{n\log_2\lambda_1+O(1)} \;\xrightarrow{n\to\infty}\; \frac{1+\log_2\lambda_1}{\log_2\lambda_1} = 1 + \frac{1}{\log_2\lambda_1}\approx 2.44,

故大 $n$ 极限为

\boxed{\;\frac{T_{\mathrm{ext}}}{T_{\mathrm{matrix}}} \;\xrightarrow{n\to\infty}\; \frac{3}{8}\cdot 2.44^{\,\alpha}.\;}

$\alpha$ 就是「乘法效率」参数： $\alpha$ 越接近 $1$ ，位数膨胀 $N'/N$ 的惩罚越轻，扩域的次数优势 $3<8$ 才显出来。

$\alpha$	乘法模型	$\frac{3}{8}\cdot 2.44^{\,\alpha}$	谁快
$1$	近线性（NTT/FFT）	$0.92$	扩域略快
$\log_2 3 \approx 1.585$	Karatsuba	$\approx 1.54$	矩阵快
$2$	朴素 $O(n^2)$	$2.23$	矩阵快

转折点由 $\frac{3}{8}\cdot 2.44^{\,\alpha} = 1$ 解出：

\alpha^\star = \frac{\log(8/3)}{\log 2.44} \approx 1.10.

$\alpha < 1.10 \rightarrow$ 扩域快； $\alpha > 1.10 \rightarrow$ 矩阵快。当然，实际的乘法模型也并非如此简单， $\alpha$ 的值也并非一直不变，这不仅是因为这个简化模型本身的局限，也是因为大整数乘法算法本身又有与 $N$ 的数量级相关的常数。

方法对比

方法	精确范围	运算模型	总代价
朴素递推	任意 $n$	$n$ 次大整数加法，第 $k$ 次位数 $\Theta(k)$	$\Theta(n^2)$ 位运算
矩阵快速幂	任意 $n$	$8\,\mu(n)$ 次 $N$ 位乘	$\Theta\bigl(8\,M(N)\,\mu(n)\bigr)$
Binet + 浮点	$n\le 75$ （本机 `f64` 实测）	$\mu(n)$ 次机器浮点乘	$\Theta(\log n)$ 浮点运算
扩域快速幂	任意 $n$	$3\,\mu(n)$ 次 $N'$ 位乘	$\Theta\bigl(3\,M(N')\,\mu(n)\bigr)$

其中 $\mu(n)=\lfloor\log_2 n\rfloor+\mathrm{popcount}(n)-1=\Theta(\log n)$ ， $N=\mathrm{bitlen}(F_n)$ ， $N'=N+n-1$ 。

浮点 Binet 在精度半径内最快；半径外用整数算法。矩阵与扩域同阶，实际快慢 $M(N)$ 与 $M(N')$ 的差距也影响很大。

本质上：矩阵法算 $A^n$ ，扩域法算 $d_1^n$ ；都是「线性递推 $\leftrightarrow$ 代数结构上的幂」。

具体实现

下面用 std.math.big.int.Managed 做精确整数，f64 做浮点 Binet。编译：

1
zig build-exe fib.zig -OReleaseFast

矩阵快速幂

1
const std = @import("std");
2
const Big = std.math.big.int.Managed;
3
const Allocator = std.mem.Allocator;
4

5
const Mat2 = struct {
6
    a: Big,
7
    b: Big,
8
    c: Big,
9
    d: Big,
10

11
    fn deinit(self: *Mat2) void {
12
        self.a.deinit();
13
        self.b.deinit();
14
        self.c.deinit();
15
        self.d.deinit();
16
    }
17

18
    fn initSet(allocator: Allocator, av: anytype, bv: anytype, cv: anytype, dv: anytype) !Mat2 {
19
        return .{
20
            .a = try Big.initSet(allocator, av),
21
            .b = try Big.initSet(allocator, bv),
22
            .c = try Big.initSet(allocator, cv),
23
            .d = try Big.initSet(allocator, dv),
24
        };
25
    }
26

27
    fn identity(allocator: Allocator) !Mat2 {
28
        return initSet(allocator, 1, 0, 0, 1);
29
    }
30

31
    fn baseA(allocator: Allocator) !Mat2 {
32
        return initSet(allocator, 1, 1, 1, 0);
33
    }
34

35
    fn copyFrom(self: *Mat2, other: *const Mat2) !void {
36
        try self.a.copy(other.a.toConst());
37
        try self.b.copy(other.b.toConst());
38
        try self.c.copy(other.c.toConst());
39
        try self.d.copy(other.d.toConst());
40
    }
41

42
    /// self = x * y，朴素 8 次整数乘
43
    fn mulInto(self: *Mat2, x: *const Mat2, y: *const Mat2, t1: *Big, t2: *Big) !void {
44
        try t1.mul(&x.a, &y.a);
45
        try t2.mul(&x.b, &y.c);
46
        try self.a.add(t1, t2);
47

48
        try t1.mul(&x.a, &y.b);
49
        try t2.mul(&x.b, &y.d);
50
        try self.b.add(t1, t2);
51

52
        try t1.mul(&x.c, &y.a);
53
        try t2.mul(&x.d, &y.c);
54
        try self.c.add(t1, t2);
55

56
        try t1.mul(&x.c, &y.b);
57
        try t2.mul(&x.d, &y.d);
58
        try self.d.add(t1, t2);
59
    }
60
};
61

62
/// F_n：A^n * (1,0)^T 的第二分量
63
fn fibMatrix(allocator: Allocator, n: u64) !Big {
64
    if (n == 0) return try Big.initSet(allocator, 0);
65
    if (n == 1) return try Big.initSet(allocator, 1);
66

67
    var result = try Mat2.identity(allocator);
68
    defer result.deinit();
69
    var base = try Mat2.baseA(allocator);
70
    defer base.deinit();
71
    var tmp = try Mat2.identity(allocator);
72
    defer tmp.deinit();
73
    var t1 = try Big.init(allocator);
74
    defer t1.deinit();
75
    var t2 = try Big.init(allocator);
76
    defer t2.deinit();
77

78
    var e = n;
79
    while (e > 0) : (e >>= 1) {
80
        if (e & 1 != 0) {
81
            try tmp.mulInto(&result, &base, &t1, &t2);
82
            try result.copyFrom(&tmp);
83
        }
84
        try tmp.mulInto(&base, &base, &t1, &t2);
85
        try base.copyFrom(&tmp);
86
    }
87
    var out = try Big.init(allocator);
88
    try out.copy(result.c.toConst());
89
    return out;
90
}

浮点 Binet

1
/// 成功时返回 F_n；超出 f64 可靠范围时返回 null（本机实测可靠到 n<=75）
2
fn fibFloat(n: u64) ?u64 {
3
    if (n == 0) return 0;
4
    if (n > 93) return null; // F_94 起 u64 装不下
5
    const lam1: f64 = (1.0 + @sqrt(5.0)) / 2.0;
6
    var base = lam1;
7
    var res: f64 = 1.0;
8
    var e = n;
9
    while (e > 0) : (e >>= 1) {
10
        if (e & 1 != 0) res *= base;
11
        base *= base;
12
    }
13
    const approx = res / @sqrt(5.0);
14
    if (!std.math.isFinite(approx)) return null;
15
    const rounded = @round(approx);
16
    if (rounded < 0.0 or rounded > @as(f64, @floatFromInt(std.math.maxInt(u64)))) return null;
17
    return @intFromFloat(rounded);
18
}

扩域 $\mathbb{Z}[\sqrt{5}]$

1
const Ext = struct {
2
    x: Big,
3
    y: Big, // 表示 x + y√5
4

5
    fn deinit(self: *Ext) void {
6
        self.x.deinit();
7
        self.y.deinit();
8
    }
9

10
    fn initSet(allocator: Allocator, xv: anytype, yv: anytype) !Ext {
11
        return .{
12
            .x = try Big.initSet(allocator, xv),
13
            .y = try Big.initSet(allocator, yv),
14
        };
15
    }
16

17
    fn copyFrom(self: *Ext, other: *const Ext) !void {
18
        try self.x.copy(other.x.toConst());
19
        try self.y.copy(other.y.toConst());
20
    }
21

22
    /// 5v = (v<<2) + v
23
    fn mul5(out: *Big, v: *const Big) !void {
24
        try out.shiftLeft(v, 2);
25
        try out.add(out, v);
26
    }
27

28
    /// self = p * q，恰 3 次大整数乘
29
    fn mulInto(self: *Ext, p: *const Ext, q: *const Ext, t1: *Big, t2: *Big, t3: *Big, t4: *Big) !void {
30
        try t1.mul(&p.x, &q.x); // ac
31
        try t2.mul(&p.y, &q.y); // bd
32
        try mul5(t4, t2); // 5 bd
33
        try self.x.add(t1, t4);
34

35
        try t3.add(&p.x, &p.y);
36
        try t4.add(&q.x, &q.y);
37
        try self.y.mul(t3, t4); // (a+b)(c+d)
38
        try t3.sub(&self.y, t1);
39
        try self.y.sub(t3, t2); // ad+bc
40
    }
41

42
    /// self = p²，恰 3 次大整数乘
43
    fn sqrInto(self: *Ext, p: *const Ext, t1: *Big, t2: *Big, t3: *Big) !void {
44
        try t1.mul(&p.x, &p.x);
45
        try t2.mul(&p.y, &p.y);
46
        try mul5(t3, t2);
47
        try self.x.add(t1, t3);
48
        try t1.mul(&p.x, &p.y);
49
        try self.y.add(t1, t1); // 2ab
50
    }
51
};
52

53
/// F_n = Y_n / 2^{n-1}
54
fn fibExt(allocator: Allocator, n: u64) !Big {
55
    if (n == 0) return try Big.initSet(allocator, 0);
56
    if (n == 1) return try Big.initSet(allocator, 1);
57

58
    var result = try Ext.initSet(allocator, 1, 0); // 1
59
    defer result.deinit();
60
    var base = try Ext.initSet(allocator, 1, 1); // 1+√5
61
    defer base.deinit();
62
    var tmp = try Ext.initSet(allocator, 0, 0);
63
    defer tmp.deinit();
64
    var t1 = try Big.init(allocator);
65
    defer t1.deinit();
66
    var t2 = try Big.init(allocator);
67
    defer t2.deinit();
68
    var t3 = try Big.init(allocator);
69
    defer t3.deinit();
70
    var t4 = try Big.init(allocator);
71
    defer t4.deinit();
72

73
    var e = n;
74
    while (e > 0) : (e >>= 1) {
75
        if (e & 1 != 0) {
76
            try tmp.mulInto(&result, &base, &t1, &t2, &t3, &t4);
77
            try result.copyFrom(&tmp);
78
        }
79
        try tmp.sqrInto(&base, &t1, &t2, &t3);
80
        try base.copyFrom(&tmp);
81
    }
82

83
    var out = try Big.init(allocator);
84
    try out.shiftRight(&result.y, @intCast(n - 1));
85
    return out;
86
}

实测

环境：Zig 0.15.2，-OReleaseFast，x86_64，CPU 为 13th Gen Intel Core i7-13700HX。计时为 wall time 平均（含调用内部分配）。

正确性：

fibFloat 与线性递推对照：精确到 $n\le 75$ ， $n=76$ 起失准；本机约 $3\text{–}5\,\mathrm{ns}$ （ $n=10\to75$ ）。
fibMatrix 与 fibExt 在 $n\in\{0,1,2,5,10,20,50,75,100,1000,10000\}$ 上结果一致。
样例： $F_{10}=55$ ， $F_{75}=2111485077978050$ ， $F_{100}=354224848179261915075$ 。

两种精确算法同为 $\Theta\bigl(M(n)\log n\bigr)$ ，区别只在常数（矩阵 $8$ 、扩域 $3$ ）与位数（ $N$ 对 $N'\approx 2.44\,N$ ）。wall-clock 几乎整段泡在大整数乘法里——两者的相对快慢由 $M$ 的增长指数 $\alpha$ 严格决定（推导见前文 §位数比的解析效应）：默认 std.math.big 走 Karatsuba（ $\alpha\approx 1.58$ ）时矩阵更快；要让「少做几次乘」变成胜势，得把 $M$ 拉近线性——下面用手写 NTT。

手写 NTT 大整数乘法

Zig 标准库没有 FFT/NTT 路径。这里用三模 NTT + CRT做卷积乘法：

素数（原根均为 $3$ ）： $998244353$ ， $469762049$ ， $167772161$
数字基 $2^{32}$ （每个 u64 limb 拆成两个 u32）

核心：蝶形与点积

1
/// 原地 Cooley–Tukey NTT。twiddles[s] 为第 s 层（len=2^{s+1}）的步长根。
2
fn nttWithTwiddles(a: []u32, mod: u32, twiddles: []const u32) void {
3
    const n = a.len;
4
    const log_n: u6 = @intCast(std.math.log2(n));
5
    bitReversePermute(a, log_n);
6

7
    var len: usize = 2;
8
    var stage: usize = 0;
9
    while (len <= n) : ({
10
        len <<= 1;
11
        stage += 1;
12
    }) {
13
        const wlen = twiddles[stage];
14
        var i: usize = 0;
15
        while (i < n) : (i += len) {
16
            var w: u32 = 1;
17
            var j: usize = 0;
18
            while (j < len / 2) : (j += 1) {
19
                const u = a[i + j];
20
                const v = modMul(a[i + j + len / 2], w, mod);
21
                a[i + j] = modAdd(u, v, mod);
22
                a[i + j + len / 2] = modSub(u, v, mod);
23
                w = modMul(w, wlen, mod);
24
            }
25
        }
26
    }
27
}
28

29
fn nttMulPoly(fa: []u32, fb: []u32, mod: u32, tw_f: []const u32, tw_i: []const u32, n_inv: u32) void {
30
    nttWithTwiddles(fa, mod, tw_f);
31
    nttWithTwiddles(fb, mod, tw_f);
32
    for (fa, fb) |*x, y| x.* = modMul(x.*, y, mod);
33
    nttWithTwiddles(fa, mod, tw_i); // 逆变换
34
    for (fa) |*x| x.* = modMul(x.*, n_inv, mod);
35
}

CRT 还原

x=\sum_{i=0}^{2} r_i\cdot\frac{M}{m_i}\cdot\Bigl(\frac{M}{m_i}\Bigr)^{-1}\bmod m_i \pmod M

1
fn crt3(r0: u32, r1: u32, r2: u32) u128 {
2
    var x: u128 = 0;
3
    const rs = [_]u32{ r0, r1, r2 };
4
    inline for (0..3) |i| {
5
        x += Mods.Mi[i] * @as(u128, rs[i]) * @as(u128, Mods.yi[i]);
6
    }
7
    return x % Mods.M;
8
}

流程：toDigits32 → 三模 nttMulPoly → 逐系数 crt3 + base- $2^{32}$ 进位 → 写回 Managed。矩阵/扩域里把 t.mul(a,b) 换成 mulNtt(t,a,b) 即可。

浮点 vs 矩阵 vs 扩域

矩阵 / 扩域底层乘法均为 NTT（mulNtt）。浮点只画精确范围 $n\le 75$ ：

浮点 vs 矩阵 vs 扩域

$n$	浮点	矩阵	扩域	扩域 / 矩阵
$10$	$3.0\,\mathrm{ns}$
$20$	$3.4\,\mathrm{ns}$
$50$	$3.9\,\mathrm{ns}$
$75$	$4.9\,\mathrm{ns}$	$1.61\,\mathrm{ms}$	$0.63\,\mathrm{ms}$	$0.39\times$
$10^{3}$	—	$2.18\,\mathrm{ms}$	$0.88\,\mathrm{ms}$	$0.40\times$
$10^{4}$	—	$5.13\,\mathrm{ms}$	$3.48\,\mathrm{ms}$	$0.68\times$
$10^{5}$	—	$39.2\,\mathrm{ms}$	$28.5\,\mathrm{ms}$	$0.73\times$
$10^{6}$	—	$443\,\mathrm{ms}$	$339\,\mathrm{ms}$	$0.76\times$

浮点是 $\Theta(\log n)$ 次机器浮点乘（纳秒级）；矩阵 / 扩域是 $\Theta(\log n)$ 次大整数乘。NTT 下扩域始终更快，但表里 扩域 / 矩阵 从 $\sim 0.4$ 抬到 $\sim 0.76$ ——相对优势在缩小。

测试代码：

递推的矩阵形式

复杂度模型

1. 矩阵快速幂

2. 特征值与对角化

特征方程

特征向量

求逆并推出通项

最近整数形式

浮点快速幂

3. 扩域

拓展：环乘法的整数乘次数

位数比的解析效应

方法对比

具体实现

矩阵快速幂

浮点 Binet

扩域 Z[5]\mathbb{Z}[\sqrt{5}]Z[5​]

实测

手写 NTT 大整数乘法

核心：蝶形与点积

CRT 还原

浮点 vs 矩阵 vs 扩域

扩域 $\mathbb{Z}[\sqrt{5}]$