参考文献:
- [Ber01] Bernstein D J. Multidigit multiplication for mathematicians[J]. Advances in Applied Mathematics, 2001: 1-19.
- FFT/NTT:以 CRT 的视角
文章目录
Map & Lift
[Ber01] 详细地调研了上百篇有关各类乘法算法的文献,以环同构的视角解释这些算法。
Mapping
映射:给定两个环 R , A R,A R,A,它们的环同态 ψ : R → A \psi: R \to A ψ:R→A 称为 “mapping R R R to A A A”
- 如果 ψ \psi ψ 是单的,那么容易求逆
- 如果 ψ \psi ψ 不是单的,一般地只需提供 r r r 的少量辅助信息,就可以从 ψ ( r ) \psi(r) ψ(r) 恢复出 r r r
例如:
-
m
∈
R
[
x
]
m \in R[x]
m∈R[x] 是首一多项式,那么
r
,
s
∈
R
[
x
]
r,s \in R[x]
r,s∈R[x] 的乘法,它可以 mapping 到商环
R
[
x
]
/
(
m
)
R[x]/(m)
R[x]/(m) 中,计算出
r
⋅
s
(
m
o
d
m
)
r \cdot s \pmod m
r⋅s(modm)
- 如果 deg m > deg r + deg s \deg m>\deg r+\deg s degm>degr+degs,那么计算结果是正确的
- 如果 deg m = deg r + deg s \deg m = \deg r+\deg s degm=degr+degs,那么只需要知道前导系数 r deg r r_{\deg r} rdegr 和 s deg s s_{\deg s} sdegs,就可以恢复出 r s = r deg r ⋅ s deg s ⋅ m + ( r ⋅ s ( m o d m ) ) rs = r_{\deg r}\cdot s_{\deg s}\cdot m + (r\cdot s \pmod m) rs=rdegr⋅sdegs⋅m+(r⋅s(modm)),这被称为 “evaluation at ∞ \infty ∞”
- 对于 r , s ∈ Z [ x ] r,s \in \mathbb Z[x] r,s∈Z[x],可以选取足够大的 M ∈ Z M \in \mathbb Z M∈Z,利用 x ↦ M x \mapsto M x↦M 环同态,计算出 r ( M ) , s ( M ) ∈ Z r(M),s(M) \in \mathbb Z r(M),s(M)∈Z,然后在整数上计算出 r ( M ) s ( M ) r(M)s(M) r(M)s(M)。根据柯西不等式,满足 M ≥ ( ∑ i r i 2 ) ( ∑ j s j 2 ) M \ge \sqrt{(\sum_i r_i^2)(\sum_j s_j^2)} M≥(∑iri2)(∑jsj2) 即可使得 r ( M ) s ( M ) = ( r s ) ( M ) r(M)s(M) = (rs)(M) r(M)s(M)=(rs)(M)
Lifting
提升:令 I I I 是环 R R R 的理想,自然的环嵌入 r + I ↦ r ′ , ∀ ( r − r ′ ) ∈ I r+I \mapsto r', \forall (r-r') \in I r+I↦r′,∀(r−r′)∈I 称为 “lifting R / I R/I R/I to R R R”
例如:
-
整数的 Base- B B B clumping:将整数环 Z \mathbb Z Z 同构映射到 Z [ y ] / ( B − y ) \mathbb Z[y]/(B-y) Z[y]/(B−y),然后提升到多项式环 Z [ y ] \mathbb Z[y] Z[y]
-
多项式的 Degree- n n n clumping(聚合):将多项式环 R [ x ] R[x] R[x] 同构映射到 R [ x ] [ y ] / ( y − x n ) R[x][y]/(y-x^n) R[x][y]/(y−xn),然后提升到 R [ x ] [ y ] R[x][y] R[x][y],它形如
f = ∑ i = 0 k n − 1 f i x i = ∑ i = 0 k − 1 ( ∑ j = 0 n − 1 f n i + j ⋅ x j ) ⋅ y i f=\sum_{i=0}^{kn-1} f_i x^i = \sum_{i=0}^{k-1}\left(\sum_{j=0}^{n-1}f_{ni+j}\cdot x^j\right) \cdot y^i f=i=0∑kn−1fixi=i=0∑k−1(j=0∑n−1fni+j⋅xj)⋅yi
将连续 n n n 个相邻系数,打包到环 R [ x ] R[x] R[x] 上的度数小于 n n n 的多项式 -
多项式的 Degree- n n n striding(跨步):将多项式环 R [ x ] R[x] R[x] 同构映射到 R [ y ] [ x ] / ( x n − y ) R[y][x]/(x^n-y) R[y][x]/(xn−y),然后提升到 R [ y ] [ x ] R[y][x] R[y][x],它形如
f = ∑ i = 0 k n − 1 f i x i = ∑ i = 0 n − 1 ( ∑ j = 0 k − 1 f i + n j ⋅ y j ) ⋅ x i f=\sum_{i=0}^{kn-1} f_i x^i = \sum_{i=0}^{n-1}\left(\sum_{j=0}^{k-1}f_{i+nj}\cdot y^j\right) \cdot x^i f=i=0∑kn−1fixi=i=0∑n−1(j=0∑k−1fi+nj⋅yj)⋅xi
将间隔 n n n 的系数,打包到环 R [ y ] R[y] R[y] 上的度数小于 k k k 的多项式
Karatsuba’s trick
为了计算
R
[
x
]
R[x]
R[x] 内的线性函数
a
=
a
0
+
a
1
x
a=a_0+a_1x
a=a0+a1x 和
b
=
b
0
+
b
1
x
b=b_0+b_1x
b=b0+b1x 的乘积,
ψ
:
R
[
x
]
→
R
[
x
]
/
(
x
2
−
x
)
≅
R
[
x
]
/
(
x
)
×
R
[
x
]
/
(
x
−
1
)
\psi: R[x] \to R[x]/(x^2-x) \cong R[x]/(x) \times R[x]/(x-1)
ψ:R[x]→R[x]/(x2−x)≅R[x]/(x)×R[x]/(x−1)
其中的 CRT 基是
{
1
−
x
,
x
}
\{1-x,x\}
{1−x,x},辅助信息
u
=
a
1
b
1
u=a_1b_1
u=a1b1
那么就有:
a
↦
(
a
0
,
a
0
+
a
1
)
b
↦
(
b
0
,
b
0
+
b
a
1
)
a
b
=
a
0
b
0
⋅
(
1
−
x
)
+
(
a
0
+
a
1
)
(
b
0
+
b
1
)
⋅
x
+
u
⋅
(
x
2
−
x
)
\begin{aligned} a &\mapsto (a_0, a_0+a_1)\\ b &\mapsto (b_0, b_0+ba_1)\\ ab &= a_0b_0 \cdot(1-x) + (a_0+a_1)(b_0+b_1) \cdot x + u \cdot (x^2-x) \end{aligned}
abab↦(a0,a0+a1)↦(b0,b0+ba1)=a0b0⋅(1−x)+(a0+a1)(b0+b1)⋅x+u⋅(x2−x)
Karatsuba multiplication for integers:给定两个整数
0
≤
a
,
b
<
2
2
n
0\le a,b < 2^{2n}
0≤a,b<22n,
- 先将 Z \mathbb Z Z 映射到 Z [ y ] / ( 2 n − y ) \mathbb Z[y]/(2^n-y) Z[y]/(2n−y),再提升到 Z [ y ] \mathbb Z[y] Z[y],此时 a ( y ) , b ( y ) a(y),b(y) a(y),b(y) 是线性多项式
- 利用 Karatsuba’s trick,计算它们的乘积 a ( y ) ⋅ b ( y ) ∈ Z [ y ] a(y) \cdot b(y) \in \mathbb Z[y] a(y)⋅b(y)∈Z[y]
- 最后带入 y = 2 n y=2^n y=2n,回到 a b ∈ Z ab \in \mathbb Z ab∈Z
Karatsuba multiplication for polynomials:给定两个多项式 a , b ∈ R [ x ] a,b \in R[x] a,b∈R[x],满足 deg a , deg b < 2 n \deg a,\deg b < 2n dega,degb<2n
- 先将 R [ x ] R[x] R[x] 度 n n n 聚合为 R [ x ] [ y ] / ( y − x n ) R[x][y]/(y-x^n) R[x][y]/(y−xn) 并提升,此时 a ( y ) , b ( y ) a(y),b(y) a(y),b(y) 是关于 y y y 线性多项式
- 利用 Karatsuba’s trick,计算它们的乘积 a ( y ) ⋅ b ( y ) ∈ ( R [ x ] ) [ y ] a(y) \cdot b(y) \in (R[x])[y] a(y)⋅b(y)∈(R[x])[y]
- 最后带入 y = x n y=x^n y=xn,回到 a b ∈ R [ x ] ab \in R[x] ab∈R[x]
Dual Karatsuba multiplication:
- 先将 R [ x ] R[x] R[x] 度 2 2 2 跨步为 R [ y ] [ x ] / ( x 2 − y ) R[y][x]/(x^2-y) R[y][x]/(x2−y) 并提升,此时 a ( x ) , b ( x ) a(x),b(x) a(x),b(x) 是关于 x x x 的线性多项式
- 利用 Karatsuba’s trick,计算它们的乘积 a ( x ) ⋅ b ( x ) ∈ ( R [ y ] ) [ x ] a(x) \cdot b(x) \in (R[y])[x] a(x)⋅b(x)∈(R[y])[x]
- 最后带入 y = x 2 y=x^2 y=x2,回到 a b ∈ R [ x ] ab \in R[x] ab∈R[x]
对于 CRT 分解后的小环 R [ x ] / ( x ) , R [ x ] / ( x − 1 ) R[x]/(x), R[x]/(x-1) R[x]/(x),R[x]/(x−1) 上的乘法,递归地使用 Karatsuba multiplication,可以达到 O ( n log 3 ) O(n^{\log3}) O(nlog3) 复杂度。这是第一个优于 n 2 n^2 n2 的算法。
Toom’s trick
Toom 推广了 Karatsuba 技巧。为了计算
R
[
x
]
R[x]
R[x] 内次数小于
k
k
k 的多项式的乘积,假设
2
,
3
,
⋯
,
2
(
k
−
1
)
2,3,\cdots,2(k-1)
2,3,⋯,2(k−1) 是环
R
R
R 上可消去的,那么
ψ
:
R
[
x
]
→
R
[
x
]
/
(
x
−
(
−
k
+
1
)
)
(
x
−
(
−
k
+
2
)
)
⋯
(
x
−
(
k
−
1
)
)
\psi: R[x] \to R[x]/(x-(-k+1))(x-(-k+2)) \cdots (x-(k-1))
ψ:R[x]→R[x]/(x−(−k+1))(x−(−k+2))⋯(x−(k−1))
继续做 CRT 分解,可以得到
2
k
−
1
2k-1
2k−1 个小环。在 CRT 域下计算乘法,然后用 CRT basis 合成,然后提升回
R
[
x
]
R[x]
R[x]
Toom multiplication for integers:对于
0
≤
a
,
b
≤
2
k
n
−
1
0\le a,b \le 2^{kn}-1
0≤a,b≤2kn−1,先执行基
2
n
2^n
2n 聚合,得到次数小于
k
k
k 的多项式,然后利用 Toom’s trick 计算多项式乘积。以
k
=
3
k=3
k=3 为例,
Z
→
Z
[
y
]
/
(
2
n
−
y
)
→
Z
[
y
]
/
(
y
5
−
5
y
3
+
4
y
)
\mathbb Z \to \mathbb Z[y]/(2^n-y) \to Z[y]/(y^5-5y^3+4y)
Z→Z[y]/(2n−y)→Z[y]/(y5−5y3+4y)
然后做 CRT 分解
Z
[
y
]
/
(
y
5
−
5
y
3
+
4
y
)
≅
∏
t
=
−
2
2
Z
[
y
]
/
(
y
−
t
)
Z[y]/(y^5-5y^3+4y) \cong \prod_{t=-2}^{2} \mathbb Z[y]/(y-t)
Z[y]/(y5−5y3+4y)≅t=−2∏2Z[y]/(y−t)
这里的小环
Z
[
y
]
/
(
y
−
t
)
≅
Z
[
2
n
]
/
(
2
n
−
t
)
\mathbb Z[y]/(y-t) \cong \mathbb Z[2^n]/(2^n-t)
Z[y]/(y−t)≅Z[2n]/(2n−t),既可以按照多项式乘法,也可以按照整数乘法。
Toom multiplication for polynomials:对于 deg a , deg b ≤ k n − 1 \deg a,\deg b \le kn-1 dega,degb≤kn−1,先执行度 n n n 聚合,得到次数小于 k k k 的多项式,然后利用 Toom’s trick 计算多项式乘积
Dual Toom multiplication:对于 deg a , deg b ≤ k n − 1 \deg a,\deg b \le kn-1 dega,degb≤kn−1,先执行度 k k k 跨步,得到次数小于 k k k 的多项式,然后利用 Toom’s trick 计算多项式乘积
对于 CRT 分解后的小环 R [ x ] / ( x − t ) , t = − k + 1 , ⋯ , k − 1 R[x]/(x-t), t=-k+1,\cdots,k-1 R[x]/(x−t),t=−k+1,⋯,k−1 上的乘法,递归地使用 Toom multiplication,可以达到 n exp ( O ( log n ) ) n\exp\left(O(\sqrt{\log n})\right) nexp(O(logn)) 复杂度。这是第一个达到 O ~ ( n ) \tilde O(n) O~(n) 的算法。
FFT trick
如果
b
∈
R
b \in R
b∈R 满足
2
b
2b
2b 是可逆的,那么下面的映射是 CRT 可逆的,
R
[
x
]
/
(
x
2
n
−
b
2
)
→
R
[
x
]
/
(
x
n
−
b
)
×
R
[
x
]
/
(
x
n
+
b
)
R[x]/(x^{2n}-b^2) \to R[x]/(x^{n}-b) \times R[x]/(x^{n}+b)
R[x]/(x2n−b2)→R[x]/(xn−b)×R[x]/(xn+b)
FFT:假如存在
ζ
∈
R
\zeta \in R
ζ∈R 满足
ζ
n
=
−
1
\zeta^{n}=-1
ζn=−1,那么递归地对
R
[
x
]
/
(
x
2
n
−
1
)
R[x]/(x^{2n}-1)
R[x]/(x2n−1) 使用 FFT trick,可以达到
O
(
n
log
n
)
O(n\log n)
O(nlogn) 复杂度。注意,只要任意的环
R
R
R 中存在
ζ
\zeta
ζ 即可,不需要它成为域。
twisted FFT:将环
R
[
x
]
/
(
x
2
n
−
1
)
R[x]/(x^{2n}-1)
R[x]/(x2n−1) 分解为
R
[
x
]
/
(
x
n
−
1
)
R[x]/(x^{n}-1)
R[x]/(xn−1) 和
R
[
x
]
/
(
x
n
+
1
)
R[x]/(x^{n}+1)
R[x]/(xn+1),然后采用同构映射
x
↦
ζ
y
x \mapsto \zeta y
x↦ζy,将后者扭曲为
R
[
x
]
/
(
x
n
+
1
)
≅
R
[
y
]
/
(
y
n
−
1
)
R[x]/(x^{n}+1) \cong R[y]/(y^{n}-1)
R[x]/(xn+1)≅R[y]/(yn−1)
接下来的递归 FFT 只需考虑形如
(
x
n
−
1
)
,
(
y
n
−
1
)
(x^n-1), (y^n-1)
(xn−1),(yn−1) 的分解,实现更加简单
split-radix FFT:假如存在
ζ
∈
R
\zeta \in R
ζ∈R 满足
ζ
2
n
=
−
1
\zeta^{2n}=-1
ζ2n=−1,简记
i
=
ζ
n
i=\zeta^n
i=ζn,执行如下形式的分解
R
[
x
]
/
(
x
4
n
−
1
)
→
R
[
x
]
/
(
x
2
n
−
1
)
×
R
[
x
]
/
(
x
n
−
i
)
×
R
[
x
]
/
(
x
n
+
i
)
R[x]/(x^{4n}-1) \to R[x]/(x^{2n}-1) \times R[x]/(x^{n}-i) \times R[x]/(x^{n}+i)
R[x]/(x4n−1)→R[x]/(x2n−1)×R[x]/(xn−i)×R[x]/(xn+i)
接着,
- 对于 R [ x ] / ( x n − i ) R[x]/(x^{n}-i) R[x]/(xn−i) 采取同构 x ↦ ζ y x \mapsto \zeta y x↦ζy,扭曲为 R [ y ] / ( y n − 1 ) R[y]/(y^n-1) R[y]/(yn−1)
- 对于 R [ x ] / ( x n + i ) R[x]/(x^{n}+i) R[x]/(xn+i) 采取同构 x ↦ ζ 3 z x \mapsto \zeta^3 z x↦ζ3z,扭曲为 R [ z ] / ( z n − 1 ) R[z]/(z^n-1) R[z]/(zn−1)
现在我们得到了
R
[
x
]
/
(
x
4
n
−
1
)
→
R
[
x
]
/
(
x
2
n
−
1
)
×
R
[
y
]
/
(
y
n
−
1
)
×
R
[
z
]
/
(
z
n
−
1
)
R[x]/(x^{4n}-1) \to R[x]/(x^{2n}-1) \times R[y]/(y^{n}-1) \times R[z]/(z^{n}-1)
R[x]/(x4n−1)→R[x]/(x2n−1)×R[y]/(yn−1)×R[z]/(zn−1)
radix-3 FFT:如果存在
w
∈
R
w \in R
w∈R 满足
1
+
w
+
w
2
=
0
1+w+w^2=0
1+w+w2=0,并且
b
∈
R
b \in R
b∈R 满足
3
b
2
3b^2
3b2 是可逆的,那么如下的分解是 CRT 可逆的,
R
[
x
]
/
(
x
3
n
−
b
3
)
→
R
[
x
]
/
(
x
n
−
b
)
×
R
[
x
]
/
(
x
n
−
w
b
)
×
R
[
x
]
/
(
x
n
−
w
2
b
)
R[x]/(x^{3n}-b^3) \to R[x]/(x^n-b) \times R[x]/(x^n-wb) \times R[x]/(x^n-w^2b)
R[x]/(x3n−b3)→R[x]/(xn−b)×R[x]/(xn−wb)×R[x]/(xn−w2b)
使用场景:
- 使用 FFT 实现 InvFFT 不是个好主意,数据移动明显偏多
- 递归形式的 FFT,各个小环的存储大小更加适配 Cache,因此一般会比迭代形式的 FFT 更快
- 特征 2 2 2 的环 R R R,应当使用 radix-3 FFT
Good’s trick
假设 gcd ( m , n ) = 1 \gcd(m,n)=1 gcd(m,n)=1,那么 Z m n ∗ ≅ Z m ∗ × Z n ∗ \mathbb Z_{mn}^* \cong \mathbb Z_m^* \times \mathbb Z_n^* Zmn∗≅Zm∗×Zn∗,映射 a → ( b , c ) a \to (b,c) a→(b,c) 的逆是 a = b ′ c ′ a=b'c' a=b′c′,其中的 b ′ , c ′ ∈ Z m n ∗ b',c' \in \mathbb Z_{mn}^* b′,c′∈Zmn∗ 是满足 b ′ ≡ b ( m o d m ) b'\equiv b \pmod m b′≡b(modm) 和 c ′ ≡ c ( m o d n ) c'\equiv c \pmod n c′≡c(modn) 的恰当代表元。
假如
x
x
x 是
m
n
mn
mn 次本原单位根,
y
,
z
y,z
y,z 分别是
m
,
n
m,n
m,n 次本原单位根,因此
y
z
yz
yz 也是一个
m
n
mn
mn 次本原单位根。映射
x
↦
y
z
x \mapsto yz
x↦yz,导致了如下的同构:
R
[
x
]
/
(
x
m
n
−
1
)
≅
(
R
[
y
]
/
(
y
m
−
1
)
)
[
z
]
/
(
z
n
−
1
)
R[x]/(x^{mn}-1) \cong (R[y]/(y^m-1))[z]/(z^n-1)
R[x]/(xmn−1)≅(R[y]/(ym−1))[z]/(zn−1)
对于
f
=
∑
i
f
i
x
i
∈
R
[
x
]
/
(
x
m
n
−
1
)
f = \sum_i f_i x^i \in R[x]/(x^{mn}-1)
f=∑ifixi∈R[x]/(xmn−1),具体的转化为:
f
(
x
)
↦
∑
i
=
0
m
n
−
1
f
i
⋅
(
y
z
)
i
=
∑
i
=
0
m
n
−
1
f
i
⋅
y
i
(
m
o
d
m
)
z
i
(
m
o
d
n
)
=
∑
i
=
0
n
−
1
(
∑
j
=
0
m
−
1
f
i
+
n
j
⋅
y
i
+
n
j
(
m
o
d
m
)
)
⋅
z
i
\begin{aligned} f(x) &\mapsto \sum_{i=0}^{mn-1} f_i \cdot (yz)^i\\ &= \sum_{i=0}^{mn-1} f_i \cdot y^{i \pmod m} z^{i \pmod n}\\ &= \sum_{i=0}^{n-1}\left(\sum_{j=0}^{m-1} f_{i+nj} \cdot y^{i+nj \pmod m}\right) \cdot z^{i}\\ \end{aligned}
f(x)↦i=0∑mn−1fi⋅(yz)i=i=0∑mn−1fi⋅yi(modm)zi(modn)=i=0∑n−1(j=0∑m−1fi+nj⋅yi+nj(modm))⋅zi
Good’s trick 往往和 FFT 组合使用,减少 Padding 长度(FFT 要求长度是二的幂次)
Manufacturing Roots of Unit
为了对 f ∈ R [ x ] f \in R[x] f∈R[x] 使用 FFT,必须要求环 R R R 中存在恰当的本原单位根。如果不存在,最简单的思路就是通过环扩张,使得在 R [ y ] R[y] R[y] 中存在,然后将系数从 R R R 提升到 R [ y ] R[y] R[y] 中,得到 f ∈ ( R [ y ] ) [ x ] f \in (R[y])[x] f∈(R[y])[x],虽然 f i ∈ R ⊆ R [ y ] f_i \in R \subseteq R[y] fi∈R⊆R[y]。
注意到如果存在 ζ ∈ R \zeta \in R ζ∈R,它已经是 n n n 次单位根,那么扩环 R [ y ] / ( y t − ζ ) R[y]/(y^t-\zeta) R[y]/(yt−ζ) 中的 y y y 就是 n t nt nt 次单位根。这通过已有的单位根,减少了环的扩张次数(如果简单使用 ( y n t − 1 ) (y^{nt}-1) (ynt−1),这扩张了 n t nt nt 次)。对于 Z [ x ] \mathbb Z[x] Z[x] 来说,系数的环可以是 Z / ( t ) \mathbb Z/(t) Z/(t)(充分大的模数 t t t,不必是素的),也可以是 C \mathbb C C(环 Z \mathbb Z Z 的扩张,但是精度受限),只要能够找出单位根,并且要求数据不发生溢出。
不过,也可以通过对多项式变换,通过把一些系数打包在一起,构造出系数的环 R [ y ] R[y] R[y](直接找出了单位根),而非简单的提升。
Schonhage-Strassen trick
对于整数环
Z
/
(
2
m
n
+
1
)
\mathbb Z/(2^{mn}+1)
Z/(2mn+1),将它同构映射为
(
Z
[
y
]
/
(
y
n
+
1
)
)
/
(
2
m
−
y
)
(\mathbb Z[y]/(y^n+1))/(2^m-y)
(Z[y]/(yn+1))/(2m−y),其中的整系数不大于
2
m
2^m
2m。我们将它提升到
Z
[
y
]
/
(
y
n
+
1
)
\mathbb Z[y]/(y^n+1)
Z[y]/(yn+1),接着可以再映射到
(
Z
/
(
2
n
k
+
1
)
)
[
y
]
/
(
y
n
+
1
)
(\mathbb Z/(2^{nk}+1))[y]/(y^n+1)
(Z/(2nk+1))[y]/(yn+1),于是:
Z
/
(
2
m
n
+
1
)
≅
Z
[
y
]
/
(
y
n
+
1
,
2
m
−
y
)
→
l
i
f
t
Z
[
y
]
/
(
y
n
+
1
)
→
m
a
p
Z
2
n
k
+
1
[
y
]
/
(
y
n
+
1
)
\mathbb Z/(2^{mn}+1) \cong \mathbb Z[y]/(y^n+1, 2^m-y) \overset{\bf lift}{\to} \mathbb Z[y]/(y^n+1) \overset{\bf map}{\to} \mathbb Z_{2^{nk}+1}[y]/(y^n+1)
Z/(2mn+1)≅Z[y]/(yn+1,2m−y)→liftZ[y]/(yn+1)→mapZ2nk+1[y]/(yn+1)
预先确定参数
m
,
n
,
k
m,n,k
m,n,k 的取值,只要满足
2
n
k
>
n
⋅
2
2
m
2^{nk} > n \cdot 2^{2m}
2nk>n⋅22m,那么
Z
/
(
2
n
k
+
1
)
\mathbb Z/(2^{nk}+1)
Z/(2nk+1) 就可以正确模拟大小为
2
m
2^m
2m 的
n
n
n 个系数的(反)卷积。
特别地,在 Z / ( 2 n k + 1 ) \mathbb Z/(2^{nk}+1) Z/(2nk+1) 中的 2 k 2^k 2k,它就是 2 n 2n 2n 次本原单位根!因此,
- 可以对 Z 2 n k + 1 [ y ] / ( y n + 1 ) \mathbb Z_{2^{nk}+1}[y]/(y^n+1) Z2nk+1[y]/(yn+1) 执行 FFT 算法,计算乘法
- 将乘积提升回 Z [ y ] / ( y n + 1 ) \mathbb Z[y]/(y^n+1) Z[y]/(yn+1),它是正确计算的
- 带入 y = 2 m y=2^m y=2m 得到 Z / ( 2 n m + 1 ) \mathbb Z/(2^{nm}+1) Z/(2nm+1) 中的最终结果
Cyclic 变体:整数环 Z / ( 2 2 m n − 1 ) \mathbb Z/(2^{2mn}-1) Z/(22mn−1),同构于 ( Z [ y ] / ( y 2 n − 1 ) ) / ( 2 m − y ) (\mathbb Z[y]/(y^{2n}-1))/(2^m-y) (Z[y]/(y2n−1))/(2m−y),提升后再映射为 ( Z [ y ] / ( y 2 n − 1 ) ) / ( y 2 n − 1 ) (\mathbb Z[y]/(y^{2n}-1))/(y^{2n}-1) (Z[y]/(y2n−1))/(y2n−1)
Schonhage’s trick
对于多项式环
R
[
x
]
/
(
x
m
n
+
1
)
R[x]/(x^{mn}+1)
R[x]/(xmn+1),做如下的同构映射
R
[
x
]
/
(
x
m
n
+
1
)
≅
R
[
x
]
[
y
]
/
(
y
n
+
1
,
x
m
−
y
)
R[x]/(x^{mn}+1) \cong R[x][y]/(y^n+1, x^{m}-y)
R[x]/(xmn+1)≅R[x][y]/(yn+1,xm−y)
它关于
x
x
x 的次数小于
m
m
m,我们将它提升到
R
[
x
]
[
y
]
/
(
y
n
+
1
)
R[x][y]/(y^n+1)
R[x][y]/(yn+1),然后再做同态映射
R
[
x
]
[
y
]
/
(
y
n
+
1
)
→
(
R
[
x
]
/
(
x
n
k
+
1
)
)
[
y
]
/
(
y
n
+
1
)
R[x][y]/(y^n+1) \to (R[x]/(x^{nk}+1))[y]/(y^n+1)
R[x][y]/(yn+1)→(R[x]/(xnk+1))[y]/(yn+1)
预先确定参数
m
,
n
,
k
m,n,k
m,n,k 的取值,只要满足
n
k
>
2
m
−
2
nk > 2m-2
nk>2m−2,那么
R
[
x
]
/
(
x
n
k
+
1
)
R[x]/(x^{nk}+1)
R[x]/(xnk+1) 就可以正确模拟长度为
m
m
m 的多项式的乘积。
特别地,在 R [ x ] / ( x n k + 1 ) R[x]/(x^{nk}+1) R[x]/(xnk+1) 中的 x k x^k xk,它就是 2 n 2n 2n 次本原单位根!因此,我们对环 ( R [ x ] / ( x n k + 1 ) ) [ y ] / ( y n + 1 ) (R[x]/(x^{nk}+1))[y]/(y^n+1) (R[x]/(xnk+1))[y]/(yn+1) 做 FFT,计算乘积之后,经过反序的 lift 和 map,回到 R [ x ] / ( x m n + 1 ) R[x]/(x^{mn}+1) R[x]/(xmn+1) 上。
Cyclic 变体:多项式环 R [ x ] / ( x 2 m n − 1 ) R[x]/(x^{2mn}-1) R[x]/(x2mn−1),同构于 R [ x ] [ y ] / ( y 2 n − 1 , x m − y ) R[x][y]/(y^{2n}-1, x^{m}-y) R[x][y]/(y2n−1,xm−y),提升后再映射为 ( R [ x ] / ( x n k + 1 ) ) [ y ] / ( y 2 n − 1 ) (R[x]/(x^{nk}+1))[y]/(y^{2n}-1) (R[x]/(xnk+1))[y]/(y2n−1)
Radix-3 变体:对于多项式环
R
[
x
]
/
(
x
2
m
n
+
x
m
n
+
1
)
R[x]/(x^{2mn}+x^{mn}+1)
R[x]/(x2mn+xmn+1),映射为
R
[
x
]
[
y
]
/
(
y
2
n
+
y
n
+
1
,
x
m
−
y
)
→
l
i
f
t
R
[
x
]
[
y
]
/
(
y
2
n
+
y
n
+
1
)
→
m
a
p
(
R
[
x
]
/
(
z
2
n
k
+
x
n
k
+
1
)
)
[
y
]
/
(
y
2
n
+
y
n
+
1
)
R[x][y]/(y^{2n}+y^{n}+1, x^m-y) \overset{\bf lift}{\to} R[x][y]/(y^{2n}+y^{n}+1) \overset{\bf map}{\to} (R[x]/(z^{2nk}+x^{nk}+1))[y]/(y^{2n}+y^{n}+1)
R[x][y]/(y2n+yn+1,xm−y)→liftR[x][y]/(y2n+yn+1)→map(R[x]/(z2nk+xnk+1))[y]/(y2n+yn+1)
其中
x
n
k
∈
R
[
x
]
/
(
z
2
n
k
+
x
n
k
+
1
)
x^{nk} \in R[x]/(z^{2nk}+x^{nk}+1)
xnk∈R[x]/(z2nk+xnk+1) 是
3
3
3 次单位根,于是就有分解
y
2
n
+
y
n
+
1
=
(
y
n
−
x
n
k
)
(
y
n
−
x
2
n
k
)
y^{2n}+y^n+1 = (y^n-x^{nk})(y^n-x^{2nk})
y2n+yn+1=(yn−xnk)(yn−x2nk)
另外
x
k
x^k
xk 是
2
n
2n
2n 次本原单位根,从而可以对它们继续执行 FFT 算法。
Nussbaumer’s trick
对于多项式环
R
[
x
]
/
(
x
m
n
k
+
1
)
R[x]/(x^{mnk}+1)
R[x]/(xmnk+1) 做跨步,得到同构
R
[
y
]
[
x
]
/
(
y
n
k
+
1
,
x
m
−
y
)
≅
(
R
[
y
]
/
(
y
n
k
+
1
)
)
[
x
]
/
(
x
m
−
y
)
R[y][x]/(y^{nk}+1, x^{m}-y) \cong (R[y]/(y^{nk}+1))[x]/(x^{m}-y)
R[y][x]/(ynk+1,xm−y)≅(R[y]/(ynk+1))[x]/(xm−y)
系数的范围是
R
[
y
]
/
(
y
n
k
+
1
)
R[y]/(y^{nk}+1)
R[y]/(ynk+1),而关于
x
x
x 的度数小于
m
m
m,将它提升后,再同态映射:
(
R
[
y
]
/
(
y
n
k
+
1
)
)
[
x
]
/
(
x
m
−
y
)
→
l
i
f
t
(
R
[
y
]
/
(
y
n
k
+
1
)
)
[
x
]
→
m
a
p
(
R
[
y
]
/
(
y
n
k
+
1
)
)
[
x
]
/
(
x
2
n
−
1
)
(R[y]/(y^{nk}+1))[x]/(x^{m}-y) \overset{\bf lift}{\to} (R[y]/(y^{nk}+1))[x] \overset{\bf map}{\to} (R[y]/(y^{nk}+1))[x]/(x^{2n}-1)
(R[y]/(ynk+1))[x]/(xm−y)→lift(R[y]/(ynk+1))[x]→map(R[y]/(ynk+1))[x]/(x2n−1)
选取参数
m
,
n
,
k
m,n,k
m,n,k,只要满足
n
≥
m
n \ge m
n≥m,那么后者就可以正确模拟
(
R
[
y
]
/
(
y
n
k
+
1
)
)
[
x
]
/
(
x
m
−
y
)
(R[y]/(y^{nk}+1))[x]/(x^{m}-y)
(R[y]/(ynk+1))[x]/(xm−y),从而计算出
R
[
x
]
/
(
x
m
n
k
+
1
)
R[x]/(x^{mnk}+1)
R[x]/(xmnk+1)
特别地, y k ∈ R [ y ] / ( y n k + 1 ) y^k \in R[y]/(y^{nk}+1) yk∈R[y]/(ynk+1) 就是 2 n 2n 2n 次本原单位根,从而可以执行 FFT 算法。
Radix-3 变体:对于多项式环
R
[
x
]
/
(
x
2
m
n
k
+
x
m
n
k
+
1
)
R[x]/(x^{2mnk}+x^{mnk}+1)
R[x]/(x2mnk+xmnk+1),映射为
(
R
[
y
]
/
(
y
2
n
k
+
y
n
k
+
1
)
)
[
x
]
/
(
x
m
−
y
)
→
l
i
f
t
(
R
[
y
]
/
(
y
2
n
k
+
y
n
k
+
1
)
)
[
x
]
→
m
a
p
(
R
[
y
]
/
(
y
2
n
k
+
y
n
k
+
1
)
)
[
x
]
/
(
x
3
n
−
1
)
(R[y]/(y^{2nk}+y^{nk}+1))[x]/(x^m-y) \overset{\bf lift}{\to} (R[y]/(y^{2nk}+y^{nk}+1))[x] \overset{\bf map}{\to} (R[y]/(y^{2nk}+y^{nk}+1))[x]/(x^{3n}-1)
(R[y]/(y2nk+ynk+1))[x]/(xm−y)→lift(R[y]/(y2nk+ynk+1))[x]→map(R[y]/(y2nk+ynk+1))[x]/(x3n−1)
其中的
y
n
k
∈
R
[
y
]
/
(
y
2
n
k
+
y
n
k
+
1
)
y^{nk} \in R[y]/(y^{2nk}+y^{nk}+1)
ynk∈R[y]/(y2nk+ynk+1) 是
3
3
3 次单位根,并且
y
k
y^{k}
yk 是
3
n
3n
3n 次单位根。
Cantor-Kaltofen theorem
对于任意的环 R R R,为了快速计算 a ⋅ b ∈ R a \cdot b \in R a⋅b∈R,可以通过 radix-2 和 radix-3 分别计算两遍,得到 a ⋅ b a \cdot b a⋅b 的缩放结果,然后由于 gcd ( 2 , 3 ) = 1 \gcd(2,3)=1 gcd(2,3)=1,从而可以将它们组合出最终结果。
令整数 e ≥ 1 , m ≥ 1 e \ge 1,m\ge 1 e≥1,m≥1,满足 2 e − 1 < m ≤ 2 e 2^{e-1}< m \le 2^e 2e−1<m≤2e,给定两个多项式 a , b ∈ R [ x ] a,b \in R[x] a,b∈R[x],那么,
第一步:令 k = ⌈ m / 2 ⌉ k=\lceil m/2 \rceil k=⌈m/2⌉,计算 K = 2 k , M = 2 m , E = 2 e K=2^k, M=2^m, E=2^e K=2k,M=2m,E=2e,定义 A = R [ y ] / ( y K + 1 ) A=R[y]/(y^K+1) A=R[y]/(yK+1)。如果 m m m 是偶数(此时 M / K = K = 2 m / 2 M/K=K=2^{m/2} M/K=K=2m/2),则 y y y 就是 2 M / K 2M/K 2M/K 次本原单位根;如果 m m m 是奇数(此时 M / K = 2 ( m − 1 ) / 2 M/K=2^{(m-1)/2} M/K=2(m−1)/2 而 K = 2 ( m + 1 ) / 2 K=2^{(m+1)/2} K=2(m+1)/2),则 y 2 y^2 y2 就是 2 M / K 2M/K 2M/K 次本原单位根。总之,本原单位根存在。
- 利用 Nussbaumer,将 R [ x ] / ( x M + 1 ) R[x]/(x^{M}+1) R[x]/(xM+1) 映射到 A [ x ] / ( x M / K − y ) A[x]/(x^{M/K}-y) A[x]/(xM/K−y),然后提升并再次映射到 A [ x ] / ( x 2 M / K − 1 ) A[x]/(x^{2M/K}-1) A[x]/(x2M/K−1)
- 利用 2 M / K 2M/K 2M/K 次本原单位根, y ∈ A y \in A y∈A 或者 y 2 ∈ A y^2 \in A y2∈A,执行环 A A A 上的 FFT 算法,将 a , b a,b a,b 转换到 FFT 域
- 在 FFT 域上计算阿达玛乘法,每个系数都是 A A A 中元素。我们递归地使用本算法,导致了缩放因子 K E / 4 KE/4 KE/4(看下一步的解释)
- 逆 FFT(忽略结束时的缩放因子),逆 Nussbaumer(数据置换)。由于环 R R R 中的因子 2 M / K 2M/K 2M/K 不一定可逆,我们只能得到缩放的结果。
- 最终,计算出 2 m + e − 1 ⋅ a b ∈ R [ x ] / ( x 2 m + 1 ) 2^{m+e-1} \cdot ab \in R[x]/(x^{2^m}+1) 2m+e−1⋅ab∈R[x]/(x2m+1),总的缩放因子是 ( 2 M / K ) ( K E / 4 ) = M E / 2 = 2 m + e − 1 (2M/K)(KE/4)=ME/2=2^{m+e-1} (2M/K)(KE/4)=ME/2=2m+e−1
第二步:令 k = ⌈ m / 2 ⌉ k=\lceil m/2 \rceil k=⌈m/2⌉,计算 K = 3 k , M = 3 m , E = 3 e K=3^k, M=3^m, E=3^e K=3k,M=3m,E=3e,定义 A = R [ y ] / ( y 2 K + y K + 1 ) A=R[y]/(y^{2K}+y^K+1) A=R[y]/(y2K+yK+1)。如果 m m m 是偶数(此时 M / K = K = 3 m / 2 M/K=K=3^{m/2} M/K=K=3m/2),则 y y y 就是 3 M / K 3M/K 3M/K 次本原单位根;如果 m m m 是奇数(此时 M / K = 3 ( m − 1 ) / 2 M/K=3^{(m-1)/2} M/K=3(m−1)/2 而 K = 3 ( m + 1 ) / 2 K=3^{(m+1)/2} K=3(m+1)/2),则 y 3 y^3 y3 就是 3 M / K 3M/K 3M/K 次本原单位根。总之,本原单位根存在。
- 利用 radix-3 Nussbaumer,将 R [ x ] / ( x 2 M + x M + 1 ) R[x]/(x^{2M}+x^M+1) R[x]/(x2M+xM+1) 映射到 A [ x ] / ( x M / K − y ) A[x]/(x^{M/K}-y) A[x]/(xM/K−y),然后提升并再次映射到 A [ x ] / ( x 2 M / K + x M / K + 1 ) A[x]/(x^{2M/K}+x^{M/K}+1) A[x]/(x2M/K+xM/K+1)
- 利用 3 M / K 3M/K 3M/K 次本原单位根, y ∈ A y \in A y∈A 或者 y 3 ∈ A y^3 \in A y3∈A,执行环 A A A 上的 radix-3 FFT 算法,将 a , b a,b a,b 转换到 FFT 域
- 在 FFT 域上计算阿达玛乘法,每个系数都是 A A A 中元素。我们递归地使用本算法,导致了缩放因子 K E / 9 KE/9 KE/9(看下一步的解释)
- 逆 FFT(忽略结束时的缩放因子),逆 Nussbaumer(数据置换)。由于环 R R R 中的因子 3 M / K 3M/K 3M/K 不一定可逆,我们只能得到缩放的结果。
- 最终,计算出 3 m + e − 1 ⋅ a b ∈ R [ x ] / ( x 3 m + 1 ) 3^{m+e-1} \cdot ab \in R[x]/(x^{3^m}+1) 3m+e−1⋅ab∈R[x]/(x3m+1),总的缩放因子是 ( 3 M / K ) ( K E / 9 ) = M E / 3 = 3 m + e − 1 (3M/K)(KE/9)=ME/3=3^{m+e-1} (3M/K)(KE/9)=ME/3=3m+e−1
第三步:给定 a , b ∈ R [ x ] a,b \in R[x] a,b∈R[x],满足 deg a + deg b < n \deg a+\deg b < n dega+degb<n,
- 我们选取参数 m , m ′ m,m' m,m′,使得两个商环的多项式度数 2 m > n , 2 × 3 m ′ > n 2^m>n, 2\times 3^{m'}>n 2m>n,2×3m′>n 足够大,计算 e = ⌈ log m ⌉ , e ′ = ⌈ log m ′ ⌉ e=\lceil \log m\rceil, e'=\lceil \log m'\rceil e=⌈logm⌉,e′=⌈logm′⌉
- 利用上述两个算法,计算出 2 m + e − 1 ⋅ a b ∈ R [ x ] / ( x 2 m + 1 ) 2^{m+e-1} \cdot ab \in R[x]/(x^{2^m}+1) 2m+e−1⋅ab∈R[x]/(x2m+1) 以及 3 m ′ + e ′ − 1 ⋅ a b ∈ R [ x ] / ( x 2 m ′ + 1 ) 3^{m'+e'-1} \cdot ab \in R[x]/(x^{2^{m'}}+1) 3m′+e′−1⋅ab∈R[x]/(x2m′+1),提升到 R [ x ] R[x] R[x]
- 因为 gcd ( 2 m + 2 − 1 , 3 m ′ + e ′ − 1 ) = 1 \gcd(2^{m+2-1},3^{m'+e'-1})=1 gcd(2m+2−1,3m′+e′−1)=1,因此存在 v < 3 m ′ + e ′ − 1 , u < 2 m + e − 1 v<3^{m'+e'-1},u<2^{m+e-1} v<3m′+e′−1,u<2m+e−1,使得 3 m ′ + e ′ − 1 u − 2 m + e − 1 v = 1 3^{m'+e'-1}u - 2^{m+e-1}v=1 3m′+e′−1u−2m+e−1v=1,从而得到 a b = u ⋅ 3 m ′ + e ′ − 1 ⋅ a b − v ⋅ 2 m + e − 1 ⋅ a b ab= u \cdot 3^{m'+e'-1} \cdot ab- v \cdot 2^{m+e-1} \cdot ab ab=u⋅3m′+e′−1⋅ab−v⋅2m+e−1⋅ab
- 本算法的复杂度为: ( 8 + 36 / log 3 ) n log n (8+36/\log3)n\log n (8+36/log3)nlogn 次基环 R R R 的乘法, ( 12 + 54 / log 3 ) n log n log ( 16 log n ) (12+54/\log3)n\log n\log(16\log n) (12+54/log3)nlognlog(16logn) 次基环 R R R 的加法,总体上是 O ( n log n ) O(n\log n) O(nlogn)
不过,虽然上述算法支持任意的环 R R R 下的 R [ x ] R[x] R[x] 中多项式的快速乘法,但是实际效率存疑。如果基环 R R R 的性质足够好(存在恰当的本原根,缩放因子可逆),那么还是 FFT 以及它的 Nussbaumer 等变体的效率更好。
Implementation
- 执行多项式环的分解时,不同 level 应当选取恰当的方法,不应局限于单一的算法
- 数据应当维持在 Transformed Versions,不要频繁执行多项式变换
- 平方的计算,相对于一般乘法,其存储开销更低。对于资源受限的,可以计算 a b = ( ( a + b ) 2 − ( a − b ) 2 ) / 4 ab=((a+b)^2-(a-b)^2)/4 ab=((a+b)2−(a−b)2)/4,用时间换空间