使用快速傅里叶变换计算大整数乘法
这是一道网络教室的题目,本来是要用分块乘法解决的,但是快速傅里叶也是一种解法。
1. 从多项式入手
实数可以表示成为多项式的形式,举一个例子
143
=
3
×
1
0
0
+
4
×
1
0
1
+
1
×
1
0
2
143=3\times10^0+4\times10^1+1\times10^2
143=3×100+4×101+1×102
这说明
143
143
143 是多项式
P
(
x
)
=
3
x
0
+
4
x
+
1
x
2
P(x)=3x^0+4x+1x^2
P(x)=3x0+4x+1x2
是当 x = 10 x=10 x=10 时的特殊情况。更一般地,所有的实数乘法都可以看作多项式相乘的形式,解决了多项式相乘的问题就解决了数乘的问题。这看似将问题复杂化了,但是不要着急,看下去就会有结果。
有两个多项式
A
(
x
)
=
∑
i
=
0
s
−
1
a
i
x
i
B
(
x
)
=
∑
i
=
0
t
−
1
b
i
x
i
A(x)=\sum_{i=0}^{s-1}a_ix^i\\ B(x)=\sum_{i=0}^{t-1}b_ix^i
A(x)=i=0∑s−1aixiB(x)=i=0∑t−1bixi
按照分配律计算
A
×
B
A\times B
A×B 则
C
(
x
)
=
∑
i
=
0
t
+
s
−
2
c
i
x
i
c
i
=
∑
j
=
0
i
∑
k
=
0
i
−
j
a
j
b
k
C(x)=\sum_{i=0}^{t+s-2}c_ix^i\\ c_i=\sum_{j=0}^{i}\sum_{k=0}^{i-j}a_jb_k
C(x)=i=0∑t+s−2cixici=j=0∑ik=0∑i−jajbk
大致的时间复杂度为 Θ ( s t ) \Theta(st) Θ(st)。
2. 换一种方式看待多项式
书接上文,多项式
A
A
A 的计算结果是由输入的变量
x
x
x 和系数
a
0
,
⋯
,
a
s
−
1
a_0,\cdots,a_{s-1}
a0,⋯,as−1 决定的,所以在表达多项式
A
A
A 时可以略去输入
x
x
x,从而看作
s
s
s 个系数。对于不同的
s
s
s 个输入
x
0
,
⋯
,
x
s
−
1
x_0,\cdots,x_{s-1}
x0,⋯,xs−1 我们可以写出
y
i
=
A
(
x
i
)
=
a
0
+
a
1
x
i
+
⋯
+
a
s
−
1
x
i
s
−
1
,
i
=
0
,
⋯
,
s
−
1
y_i=A(x_i)=a_0+a_1x_i+\cdots+a_{s-1}x^{s-1}_i,i=0,\cdots,s-1
yi=A(xi)=a0+a1xi+⋯+as−1xis−1,i=0,⋯,s−1
这是可以写成矩阵形式的
[
1
x
0
x
0
2
⋯
x
0
s
−
1
1
x
1
x
1
2
⋯
x
1
s
−
1
⋮
1
x
s
−
1
x
s
−
1
2
⋯
x
s
−
1
s
−
1
]
[
a
0
a
1
⋮
a
s
−
1
]
=
[
y
0
y
1
⋮
y
s
−
1
]
\begin{bmatrix} 1&x_0&x_0^2&\cdots&x_0^{s-1}\\ 1&x_1&x_1^2&\cdots&x_1^{s-1}\\ &&\vdots\\ 1&x_{s-1}&x_{s-1}^2&\cdots&x_{s-1}^{s-1} \end{bmatrix} \begin{bmatrix} a_0\\a_1\\\vdots\\a_{s-1} \end{bmatrix}= \begin{bmatrix} y_0\\y_1\\\vdots\\y_{s-1} \end{bmatrix}
⎣⎢⎢⎢⎡111x0x1xs−1x02x12⋮xs−12⋯⋯⋯x0s−1x1s−1xs−1s−1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮as−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y0y1⋮ys−1⎦⎥⎥⎥⎤
由于范德蒙行列式
∣
1
x
0
x
0
2
⋯
x
0
s
−
1
1
x
1
x
1
2
⋯
x
1
s
−
1
⋮
1
x
s
−
1
x
s
−
1
2
⋯
x
s
−
1
s
−
1
∣
≠
0
\begin{vmatrix} 1&x_0&x_0^2&\cdots&x_0^{s-1}\\ 1&x_1&x_1^2&\cdots&x_1^{s-1}\\ &&\vdots\\ 1&x_{s-1}&x_{s-1}^2&\cdots&x_{s-1}^{s-1} \end{vmatrix}\neq0
∣∣∣∣∣∣∣∣∣111x0x1xs−1x02x12⋮xs−12⋯⋯⋯x0s−1x1s−1xs−1s−1∣∣∣∣∣∣∣∣∣=0
故对于多项式
A
A
A 即使在不知道系数的取值但是知道
x
i
,
y
i
x_i,y_i
xi,yi 的取值的情况下也可以通过矩阵逆运算得出系数的值,
a
=
V
−
1
y
a=V^{-1}y
a=V−1y。所以多项式
A
A
A 也可以表达为自变量取值-多项式值对的形式,称为点值表达。
A
→
{
(
x
0
,
y
0
)
,
⋯
,
(
x
s
−
1
,
y
s
−
1
)
}
B
→
{
(
x
0
,
y
0
′
)
,
⋯
,
(
x
t
−
1
,
y
t
−
1
′
)
}
A\to\{(x_0,y_0),\cdots,(x_{s-1},y_{s-1})\}\\ B\to\{(x_0,y_0'),\cdots,(x_{t-1},y_{t-1}')\}
A→{(x0,y0),⋯,(xs−1,ys−1)}B→{(x0,y0′),⋯,(xt−1,yt−1′)}
这样表达带来了一个优势,那就是乘法的速度变快了。我们假设
s
=
t
=
m
s=t=m
s=t=m(即使
s
≠
t
s\neq t
s=t也可以通过添加零系数的方式添到一样),那么
C
→
{
(
x
0
,
y
0
y
0
′
)
,
⋯
,
(
x
n
−
1
,
y
m
−
1
y
m
−
1
′
)
}
C\to\{(x_0,y_0y_0'),\cdots,(x_{n-1},y_{m-1}y_{m-1}')\}
C→{(x0,y0y0′),⋯,(xn−1,ym−1ym−1′)}
的计算时间复杂度在知道 A , B A,B A,B 的点值表达的情况下为 Θ ( n ) \Theta(n) Θ(n)。但是我们需要注意对于多项式 C C C 拥有最多 2 m 2m 2m 项系数,所以如果需要通过矩阵求逆 c = V − 1 y c=V^{-1}y c=V−1y 的方式求出多项式 C C C 的常规表达法,则需要 2 m 2m 2m 个不同的 x x x 取值。
那么现在任务就明确了,对于两个输入的大整数(A, B)
- 我们按照位数写成多项式的形式,然后求出对应且足够的点值表达
- 然后计算出相乘后的点值表达(C),时间复杂度为 Θ ( m ) \Theta(m) Θ(m)
- 最后通过矩阵求逆的运算计算出点值表达对应的系数
- 对系数进行整理从而得到整数乘法的结果。
3. 单位复根
为了加速运算,我们引入单位复根这一概念。它的使用将在后面看到。
n次单位复根:使得 w n = 1 , n ∈ N w^n=1,n\in N wn=1,n∈N 的复数 w w w.
主n次单位复根: w n = exp ( 2 π i / n ) w_n=\exp(2\pi i/n) wn=exp(2πi/n).
这里给出几个简单的性质以供使用和参考:
消去定律: w d n d k = w n k . w_{dn}^{dk}=w_n^k. wdndk=wnk.
折半定律: w n n / 2 = w 2 = − 1 , n is even. w_{n}^{n/2}=w_2=-1,n\text{ is even.} wnn/2=w2=−1,n is even.
求和定律:对于任意整数
n
⩾
1
n\geqslant1
n⩾1 和非负整数
k
,
gcd
(
k
,
n
)
=
1
k,\gcd(k,n)=1
k,gcd(k,n)=1,有
∑
j
=
0
n
−
1
(
w
n
k
)
j
=
0.
\sum_{j=0}^{n-1}(w_n^k)^j=0.
j=0∑n−1(wnk)j=0.
求和定律的证明可以使用等比数列求和得到。
4. 求点值表达
对于多项式 A , B A,B A,B,我们首先通过添加零系数的方式将它们的项数补充到一致,然后再通过添加零系数的方式将两个多项式的项数都添加到一个 2 2 2 的次幂的数上,将这个数翻倍,设为 n n n。这里翻倍的原因是为了使得方程有解,原因在前文中已经说过。
令
x
i
=
w
n
i
x_i=w_n^i
xi=wni,我们要求出
A
(
w
n
0
)
,
A
(
w
n
1
)
,
⋯
,
A
(
w
n
n
−
1
)
B
(
w
n
0
)
,
B
(
w
n
1
)
,
⋯
,
B
(
w
n
n
−
1
)
A(w_n^0),A(w_n^1),\cdots,A(w_n^{n-1})\\B(w_n^0),B(w_n^1),\cdots,B(w_n^{n-1})
A(wn0),A(wn1),⋯,A(wnn−1)B(wn0),B(wn1),⋯,B(wnn−1)
以下以
A
(
x
)
A(x)
A(x) 为例,由于扩充后
A
A
A 的项数为
2
2
2 的次幂,我们设
A
[
0
]
=
a
0
+
a
2
x
+
⋯
+
a
n
−
2
x
n
/
2
−
1
A
[
1
]
=
a
1
+
a
3
x
+
⋯
+
a
n
−
1
x
n
/
2
−
1
A^{[0]}=a_0+a_2x+\cdots+a_{n-2}x^{n/2-1}\\ A^{[1]}=a_1+a_3x+\cdots+a_{n-1}x^{n/2-1}
A[0]=a0+a2x+⋯+an−2xn/2−1A[1]=a1+a3x+⋯+an−1xn/2−1
则有 A ( x ) = A [ 0 ] ( x 2 ) + x A [ 1 ] ( x 2 ) A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2) A(x)=A[0](x2)+xA[1](x2),给出其计算复杂度 T ( n ) = 2 T ( n / 2 ) + Θ ( n ) = Θ ( n lg n ) T(n)=2T(n/2)+\Theta(n)=\Theta(n\lg n) T(n)=2T(n/2)+Θ(n)=Θ(nlgn)。
特别说明一下,设
y
k
,
y
k
[
0
]
,
y
k
[
1
]
,
k
=
0
,
⋯
,
n
/
2
−
1
y_k,y_k^{[0]},y_k^{[1]},k=0,\cdots,n/2-1
yk,yk[0],yk[1],k=0,⋯,n/2−1 为多项式
A
(
w
n
k
)
,
A
[
0
]
(
w
n
2
k
)
,
A
[
1
]
(
w
n
2
k
)
A(w_n^k),A^{[0]}(w_n^{2k}),A^{[1]}(w_n^{2k})
A(wnk),A[0](wn2k),A[1](wn2k) 的值,那么
y
k
=
y
k
[
0
]
+
w
n
k
y
k
[
0
]
=
A
[
0
]
(
w
n
2
k
)
+
w
n
k
A
[
1
]
(
w
n
2
k
)
=
A
(
w
n
k
)
y
k
+
n
/
2
=
y
k
[
0
]
−
w
n
k
y
k
[
0
]
=
y
k
[
0
]
+
w
n
k
+
n
/
2
y
k
[
0
]
=
A
[
0
]
(
w
n
2
k
+
n
)
+
w
n
k
+
n
/
2
A
[
1
]
(
w
n
2
k
+
n
)
=
A
[
0
]
(
w
n
2
k
)
−
w
n
k
A
[
1
]
(
w
n
2
k
)
=
A
(
w
n
k
+
n
/
2
)
\begin{aligned} y_k&=y_k^{[0]}+w_n^ky_k^{[0]}\\ &=A^{[0]}(w_n^{2k})+w_n^kA^{[1]}(w_n^{2k})\\ &=A(w_n^k)\\ y_{k+n/2}&=y_k^{[0]}-w_n^ky_k^{[0]}\\ &=y_k^{[0]}+w_n^{k+n/2}y_k^{[0]}\\ &=A^{[0]}(w_n^{2k+n})+w_n^{k+n/2}A^{[1]}(w_n^{2k+n})\\ &=A^{[0]}(w_n^{2k})-w_n^kA^{[1]}(w_n^{2k})\\ &=A(w_n^{k+n/2}) \end{aligned}
ykyk+n/2=yk[0]+wnkyk[0]=A[0](wn2k)+wnkA[1](wn2k)=A(wnk)=yk[0]−wnkyk[0]=yk[0]+wnk+n/2yk[0]=A[0](wn2k+n)+wnk+n/2A[1](wn2k+n)=A[0](wn2k)−wnkA[1](wn2k)=A(wnk+n/2)
这指导我们我们只需遍历 k = 0 , 1 , ⋯ , n / 2 − 1 k=0,1,\cdots,n/2-1 k=0,1,⋯,n/2−1 就可以计算出 A ( w n 0 ) , ⋯ , A ( w n n − 1 ) A(w_n^0),\cdots,A(w_n^{n-1}) A(wn0),⋯,A(wnn−1)。
5. 求系数
通过上述方式我们可以求出多项式
C
=
A
×
B
C=A\times B
C=A×B 的点表达式,且时间复杂度为
Θ
(
n
lg
n
)
\Theta(n\lg n)
Θ(nlgn),接下来需要通过点表达式求出各项系数,这关键是要求出逆矩阵乘法的结果,即
[
c
0
⋮
c
n
−
1
]
=
V
(
w
n
0
,
w
n
1
,
⋯
,
w
n
n
−
1
)
−
1
[
y
0
⋮
y
n
−
1
]
\begin{bmatrix} c_0\\\vdots\\c_{n-1} \end{bmatrix}= V(w_n^0,w_n^1,\cdots,w_n^{n-1})^{-1}\begin{bmatrix}y_0\\\vdots\\y_{n-1}\end{bmatrix}
⎣⎢⎡c0⋮cn−1⎦⎥⎤=V(wn0,wn1,⋯,wnn−1)−1⎣⎢⎡y0⋮yn−1⎦⎥⎤
由于单位复根的性质,可以证明 V ( w n 0 , w n 1 , ⋯ , w n n − 1 ) − 1 = [ w n − j k / n ] n × n V(w_n^0,w_n^1,\cdots,w_n^{n-1})^{-1}=[w_n^{-jk}/n]_{n\times n} V(wn0,wn1,⋯,wnn−1)−1=[wn−jk/n]n×n,即这个范德蒙矩阵的逆矩阵的第 j j j 行第 k k k 列的元素为 w n − j k / n w^{-jk}_n/n wn−jk/n。
相当于
[
c
0
⋮
c
n
−
1
]
=
1
n
V
(
w
n
0
,
w
n
−
1
,
⋯
,
w
n
−
(
n
−
1
)
)
[
y
0
⋮
y
n
−
1
]
\begin{bmatrix} c_0\\\vdots\\c_{n-1} \end{bmatrix}= \frac1nV(w_n^0,w_n^{-1},\cdots,w_n^{-(n-1)})\begin{bmatrix}y_0\\\vdots\\y_{n-1}\end{bmatrix}
⎣⎢⎡c0⋮cn−1⎦⎥⎤=n1V(wn0,wn−1,⋯,wn−(n−1))⎣⎢⎡y0⋮yn−1⎦⎥⎤
利用之前的方法同样可以在 Θ ( n lg n ) \Theta(n\lg n) Θ(nlgn) 的时间内完成运算。
6. 例子
这里举例说明以上的运算方式
a
=
123
,
b
=
89
a
×
b
a=123,b=89\\ a\times b
a=123,b=89a×b
构建出两个多项式
A
(
x
)
=
3
+
2
x
+
x
2
,
B
=
9
+
8
x
A(x)=3+2x+x^2,B=9+8x
A(x)=3+2x+x2,B=9+8x
根据位数,最终的结果的位数有可能达到
3
+
2
=
7
3+2=7
3+2=7 位,向上填充到
2
2
2 的次幂数
8
8
8,所以我们需要8次单位复根。
exp
(
i
π
/
4
×
0
)
,
exp
(
i
π
/
4
×
1
)
,
⋯
,
exp
(
i
π
/
4
×
7
)
\exp(i\pi/4\times0),\exp(i\pi/4\times1),\cdots,\exp(i\pi/4\times7)
exp(iπ/4×0),exp(iπ/4×1),⋯,exp(iπ/4×7)
分别带入 A , B A,B A,B 得到 C C C
w | A | B | C |
---|---|---|---|
w 8 0 w_8^0 w80 | 6.000 + 0.000 i 6.000+0.000i 6.000+0.000i | 17.000 + 0.000 i 17.000+0.000i 17.000+0.000i | 102.000 + 0.000 i 102.000+0.000i 102.000+0.000i |
w 8 1 w_8^1 w81 | 4.414 + 2.414 i 4.414+2.414i 4.414+2.414i | 14.657 + 5.657 i 14.657+5.657i 14.657+5.657i | 51.042 + 60.355 i 51.042+60.355i 51.042+60.355i |
w 8 2 w_8^2 w82 | 2.000 + 2.000 i 2.000+2.000i 2.000+2.000i | 9.000 + 8.000 i 9.000+8.000i 9.000+8.000i | 2.000 + 34.000 i 2.000+34.000i 2.000+34.000i |
w 8 3 w_8^3 w83 | 1.586 + 0.414 i 1.586+0.414i 1.586+0.414i | 3.343 + 5.657 i 3.343+5.657i 3.343+5.657i | 2.958 + 10.355 i 2.958+10.355i 2.958+10.355i |
w 8 4 w_8^4 w84 | 2.000 + 0.000 i 2.000+0.000i 2.000+0.000i | 1.000 + 0.000 i 1.000+0.000i 1.000+0.000i | 2.000 + 0.000 i 2.000+0.000i 2.000+0.000i |
w 8 5 w_8^5 w85 | 1.586 − 0.414 i 1.586-0.414i 1.586−0.414i | 3.343 − 5.657 i 3.343-5.657i 3.343−5.657i | 2.958 − 10.355 i 2.958-10.355i 2.958−10.355i |
w 8 6 w_8^6 w86 | 2.000 − 2.000 i 2.000-2.000i 2.000−2.000i | 9.000 − 8.000 i 9.000-8.000i 9.000−8.000i | 2.000 − 34.000 i 2.000-34.000i 2.000−34.000i |
w 8 7 w_8^7 w87 | 4.414 − 2.414 i 4.414-2.414i 4.414−2.414i | 14.657 − 5.657 i 14.657-5.657i 14.657−5.657i | 51.042 − 60.355 i 51.042-60.355i 51.042−60.355i |
所以对于多项式
C
C
C 的点值表达有
[
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
1.000
+
0.000
i
0.707
+
0.707
i
0.000
+
1.000
i
−
0.707
+
0.707
i
−
1.000
+
0.000
i
−
0.707
−
0.707
i
−
0.000
−
1.000
i
0.707
−
0.707
i
1.000
+
0.000
i
0.000
+
1.000
i
−
1.000
+
0.000
i
−
0.000
−
1.000
i
1.000
+
0.000
i
0.000
+
1.000
i
−
1.000
+
0.000
i
−
0.000
−
1.000
i
1.000
+
0.000
i
−
0.707
+
0.707
i
−
0.000
−
1.000
i
0.707
+
0.707
i
−
1.000
+
0.000
i
0.707
−
0.707
i
0.000
+
1.000
i
−
0.707
−
0.707
i
1.000
+
0.000
i
−
1.000
+
0.000
i
1.000
+
0.000
i
−
1.000
+
0.000
i
1.000
+
0.000
i
−
1.000
+
0.000
i
1.000
+
0.000
i
−
1.000
+
0.000
i
1.000
+
0.000
i
−
0.707
−
0.707
i
0.000
+
1.000
i
0.707
−
0.707
i
−
1.000
+
0.000
i
0.707
+
0.707
i
−
0.000
−
1.000
i
−
0.707
+
0.707
i
1.000
+
0.000
i
−
0.000
−
1.000
i
−
1.000
+
0.000
i
0.000
+
1.000
i
1.000
+
0.000
i
−
0.000
−
1.000
i
−
1.000
+
0.000
i
0.000
+
1.000
i
1.000
+
0.000
i
0.707
−
0.707
i
−
0.000
−
1.000
i
−
0.707
−
0.707
i
−
1.000
+
0.000
i
−
0.707
+
0.707
i
0.000
+
1.000
i
0.707
+
0.707
i
]
[
c
0
c
1
c
2
c
3
c
4
c
5
c
6
c
7
]
=
[
102.000
+
0.000
i
51.042
+
60.355
i
2.000
+
34.000
i
2.958
+
10.355
i
2.000
+
0.000
i
2.958
−
10.355
i
2.000
−
34.000
i
51.042
−
60.355
i
]
\begin{bmatrix} 1.000+0.000i&1.000+0.000i&1.000+0.000i&1.000+0.000i&1.000+0.000i&1.000+0.000i&1.000+0.000i&1.000+0.000i\\ 1.000+0.000i&0.707+0.707i&0.000+1.000i&-0.707+0.707i&-1.000+0.000i&-0.707-0.707i&-0.000-1.000i&0.707-0.707i\\ 1.000+0.000i&0.000+1.000i&-1.000+0.000i&-0.000-1.000i&1.000+0.000i&0.000+1.000i&-1.000+0.000i&-0.000-1.000i\\ 1.000+0.000i&-0.707+0.707i&-0.000-1.000i&0.707+0.707i&-1.000+0.000i&0.707-0.707i&0.000+1.000i&-0.707-0.707i\\ 1.000+0.000i&-1.000+0.000i&1.000+0.000i&-1.000+0.000i&1.000+0.000i&-1.000+0.000i&1.000+0.000i&-1.000+0.000i\\ 1.000+0.000i&-0.707-0.707i&0.000+1.000i&0.707-0.707i&-1.000+0.000i&0.707+0.707i&-0.000-1.000i&-0.707+0.707i\\ 1.000+0.000i&-0.000-1.000i&-1.000+0.000i&0.000+1.000i&1.000+0.000i&-0.000-1.000i&-1.000+0.000i&0.000+1.000i\\ 1.000+0.000i&0.707-0.707i&-0.000-1.000i&-0.707-0.707i&-1.000+0.000i&-0.707+0.707i&0.000+1.000i&0.707+0.707i\\ \end{bmatrix} \begin{bmatrix} c_0\\ c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7 \end{bmatrix} =\begin{bmatrix} 102.000+0.000i\\ 51.042+60.355i\\ 2.000+34.000i\\ 2.958+10.355i\\ 2.000+0.000i\\ 2.958-10.355i\\ 2.000-34.000i\\ 51.042-60.355i \end{bmatrix}
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i1.000+0.000i0.707+0.707i0.000+1.000i−0.707+0.707i−1.000+0.000i−0.707−0.707i−0.000−1.000i0.707−0.707i1.000+0.000i0.000+1.000i−1.000+0.000i−0.000−1.000i1.000+0.000i0.000+1.000i−1.000+0.000i−0.000−1.000i1.000+0.000i−0.707+0.707i−0.000−1.000i0.707+0.707i−1.000+0.000i0.707−0.707i0.000+1.000i−0.707−0.707i1.000+0.000i−1.000+0.000i1.000+0.000i−1.000+0.000i1.000+0.000i−1.000+0.000i1.000+0.000i−1.000+0.000i1.000+0.000i−0.707−0.707i0.000+1.000i0.707−0.707i−1.000+0.000i0.707+0.707i−0.000−1.000i−0.707+0.707i1.000+0.000i−0.000−1.000i−1.000+0.000i0.000+1.000i1.000+0.000i−0.000−1.000i−1.000+0.000i0.000+1.000i1.000+0.000i0.707−0.707i−0.000−1.000i−0.707−0.707i−1.000+0.000i−0.707+0.707i0.000+1.000i0.707+0.707i⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡c0c1c2c3c4c5c6c7⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡102.000+0.000i51.042+60.355i2.000+34.000i2.958+10.355i2.000+0.000i2.958−10.355i2.000−34.000i51.042−60.355i⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
求解得到
[
c
0
c
1
c
2
c
3
c
4
c
5
c
6
c
7
]
=
[
27
42
25
8
0
0
0
0
]
\begin{bmatrix} c_0\\ c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7 \end{bmatrix}=\begin{bmatrix} 27\\ 42\\ 25\\ 8\\ 0\\ 0\\ 0\\ 0 \end{bmatrix}
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡c0c1c2c3c4c5c6c7⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡27422580000⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
即 C ( x ) = 27 + 42 x + 25 x 2 + 8 x 3 = A ( x ) × B ( x ) ⇒ C ( 10 ) = A ( 10 ) × B ( 10 ) ⇒ 27 + 42 × 10 + 25 × 1 0 2 + 8 × 1 0 3 = 10947 = 123 × 89 C(x)=27+42x+25x^2+8x^3=A(x)\times B(x)\Rightarrow C(10)=A(10)\times B(10)\Rightarrow 27+42\times10+25\times10^2+8\times10^3=10947=123\times 89 C(x)=27+42x+25x2+8x3=A(x)×B(x)⇒C(10)=A(10)×B(10)⇒27+42×10+25×102+8×103=10947=123×89.
7. 代码
#include <iostream>
#include <complex>
#include <cmath>
#include <vector>
using namespace std;
typedef complex<double> dcomp;
class mypoly {
public:
size_t len;
int *memory;
explicit mypoly(size_t l) : memory(new int[l]), len(l) {}
mypoly(mypoly &&x) noexcept {
len = x.len;
memory = x.memory;
x.memory = nullptr;
}
~mypoly() {
delete[] memory;
}
void read(string str) const {
for (size_t i = 0; i < str.length(); i++)
memory[i] = str[str.length() - i - 1] - '0';
for (size_t i = str.length(); i < len; i++)
memory[i] = 0;
}
};
// 正向傅里叶变换
vector<dcomp> fft(const mypoly &a, size_t off, size_t sep) {
if (sep == a.len) {
dcomp a0 = a.memory[off];
vector<dcomp> r;
r.push_back(a0);
return r;
} else {
size_t n = a.len / sep;
dcomp wn(cos(2 * M_PI / n), sin(2 * M_PI / n)), w(1., 0.);
vector<dcomp> y0 = fft(a, off, sep << 1)
, y1 = fft(a, off + sep, sep << 1)
, y(a.len / sep);
for (size_t k = 0; k < y.size() >> 1; k++) {
y[k] = y0[k] + w * y1[k];
y[k + (n >> 1)] = y0[k] - w * y1[k];
w = w * wn;
}
return y;
}
}
// 反向傅里叶变换
vector<dcomp> ifft(const vector<dcomp> &y, size_t off, size_t sep) {
if (sep == y.size()) {
dcomp y0 = y[off];
vector<dcomp> r;
r.push_back(y0);
return r;
} else {
size_t n = y.size() / sep;
dcomp wn(cos(2 * M_PI / n), -sin(2 * M_PI / n)), w(1., 0.);
vector<dcomp> a0 = ifft(y, off, sep << 1)
, a1 = ifft(y, off + sep, sep << 1)
, a(y.size() / sep);
for (size_t k = 0; k < a.size() >> 1; k++) {
a[k] = a0[k] + w * a1[k];
a[k + (n >> 1)] = a0[k] - w * a1[k];
w = w * wn;
}
return a;
}
}
int main() {
string str1, str2;
cin >> str1 >> str2;
size_t m = str1.size() | str2.size(), l = 1;
while (l <= m) l <<= 1;
l <<= 1;
mypoly p1(l), p2(l);
p1.read(str1), p2.read(str2);
vector<dcomp> r1 = fft(p1, 0, 1), r2 = fft(p2, 0, 1), r(l);
for (size_t ind = 0; ind < l; ind++)
r[ind] = r1[ind] * r2[ind];
vector<dcomp> a = ifft(r, 0, 1);
vector<int> out(a.size());
for (size_t k = 0; k < out.size(); k++)
// 因为整数截断的原因所以加0.5再向下取整防止诸如1.9998被截断成1而不是2
out[k] = (int)(std::abs(a[k]) / l + .5);
for (size_t k = 0; k < out.size() - 1; k++) {
out[k + 1] += out[k] / 10;
out[k] %= 10;
}
size_t end = out.size();
while(out[--end] == 0);
while(end != 0)
cout << out[end--];
cout << out[0] << endl;
return 0;
}
匆匆写就,如有错误恳请指出。