FFT
可以用 O(
)的時間將一個用係數表示的多項式轉換成它的點值表示的演算法
多項式有兩種表示方法:係數表示法和點值表示法。
對於多項式
係數表示法:用每一項的係數來表示
,即
點值表示法
:用平面直角座標系中n個不同的點來表示一個多項式。
證明:
那麼係數矩陣設為B=
矩陣的行列式等於矩陣轉置的行列式。將B轉置即位範德蒙行列式,因此B的行列式|B|為:
所以B滿秩,存在唯一解,所以任意不同的n個點可以唯一確定一個n次多項式。
點表示法乘法:
現在的主要問題是如何快速的將多項式的係數表示法轉為點表示法進行運算,然後快速的轉為係數表示法,DFT可以在
時間內將係數表示法轉為為點表示法,IDFT可以在
時間內將點表示法轉為係數表示法。
DFT
複數:定義
=−1,一個複數 z 可以表示為z=a+bi (a,b∈R ) 其中a稱為實部,b稱為虛部,i稱為虛數單位。
同時可以把複數看成複平面直角座標系上的一個點,如圖所示:
其中x是實數集中的座標軸,y軸是虛數單位i的座標軸。
這個點 (2,3)表示的複數就是 2+3i或者想象它代表的向量為(2,3)。
一個複數z的模定義為它到原點的距離,記為
。
複數的運算:
進入正題:
為了方便求解n一般取2的整次冪。如果暴力的找n個點帶入多項式求解,複雜度將達到
,但如果帶入一些特殊的點將大大運算。
複數域單位圓上所有的點任何次方的結果模長都為1,我們可以不用做全部的次方運算。
考慮將圓八等分:
橙色點即為n=8時要取的點,從(1,0)點開始,逆時針從0號開始標號,標到 7號
記編號為k的點代表的複數值為
,根據模長相乘,極角相加可知
,其中
稱為n次單位根,每個
都可以求出:
用圖表示為:
即為要代入的
。
單位根的性質:
1。
證明:
2。
證明:
3。
簡化運算證明:
觀察上述式子可知:如果已知
和
可以同時知道
和
。因此可以用遞迴在
時間內完成DFT運算。
程式碼如圖所示:
IDFT
現在快速把點值表示式轉化為係數表示式。
於是我們發現只需要把
換成
,然後DFT,然後將所得的結果除以n,就可以實現IDFT了。
遞迴版FFT程式碼:
void
fft
(
int
n
,
complex
<
double
>*
buffer
,
int
offset
,
int
step
,
complex
<
double
>*
epsilon
)
{
if
(
n
==
1
)
return
;
int
m
=
n
>>
1
;
fft
(
m
,
buffer
,
offset
,
step
<<
1
,
epsilon
);
fft
(
m
,
buffer
,
offset
+
step
,
step
<<
1
,
epsilon
);
for
(
int
k
=
0
;
k
!=
m
;
++
k
)
{
int
pos
=
2
*
step
*
k
;
temp
[
k
]
=
buffer
[
pos
+
offset
]
+
epsilon
[
k
*
step
]
*
buffer
[
pos
+
offset
+
step
];
temp
[
k
+
m
]
=
buffer
[
pos
+
offset
]
-
epsilon
[
k
*
step
]
*
buffer
[
pos
+
offset
+
step
];
}
for
(
int
i
=
0
;
i
!=
n
;
++
i
)
buffer
[
i
*
step
+
offset
]
=
temp
[
i
];
}
迭代版FFT和蝴蝶操作:
遞迴可以實現FFT,但常數過大,考慮迭代最佳化操作。
接下來考慮如果翻轉二進位制數:
#include
#include
#include
#include
using
namespace
std
;
const
int
N
=
300010
;
const
double
PI
=
acos
(
-
1
);
int
n
,
m
;
struct
Complex
{
double
x
,
y
;
Complex
operator
+
(
const
Complex
&
t
)
const
{
return
{
x
+
t
。
x
,
y
+
t
。
y
};
}
Complex
operator
-
(
const
Complex
&
t
)
const
{
return
{
x
-
t
。
x
,
y
-
t
。
y
};
}
Complex
operator
*
(
const
Complex
&
t
)
const
{
return
{
x
*
t
。
x
-
y
*
t
。
y
,
x
*
t
。
y
+
y
*
t
。
x
};
}
}
a
[
N
],
b
[
N
];
int
rev
[
N
],
bit
,
tot
;
void
fft
(
Complex
a
[],
int
inv
)
{
for
(
int
i
=
0
;
i
<
tot
;
i
++
)
if
(
i
<
rev
[
i
])
swap
(
a
[
i
],
a
[
rev
[
i
]]);
for
(
int
mid
=
1
;
mid
<
tot
;
mid
<<=
1
)
//從下到上迭代
{
auto
w1
=
Complex
({
cos
(
PI
/
mid
),
inv
*
sin
(
PI
/
mid
)});
for
(
int
i
=
0
;
i
<
tot
;
i
+=
mid
*
2
)
{
auto
wk
=
Complex
({
1
,
0
});
for
(
int
j
=
0
;
j
<
mid
;
j
++
,
wk
=
wk
*
w1
)
{
auto
x
=
a
[
i
+
j
],
y
=
wk
*
a
[
i
+
j
+
mid
];
a
[
i
+
j
]
=
x
+
y
,
a
[
i
+
j
+
mid
]
=
x
-
y
;
}
}
}
}
int
main
()
{
scanf
(
“%d%d”
,
&
n
,
&
m
);
for
(
int
i
=
0
;
i
<=
n
;
i
++
)
scanf
(
“%lf”
,
&
a
[
i
]。
x
);
for
(
int
i
=
0
;
i
<=
m
;
i
++
)
scanf
(
“%lf”
,
&
b
[
i
]。
x
);
while
((
1
<<
bit
)
<
n
+
m
+
1
)
bit
++
;
tot
=
1
<<
bit
;
for
(
int
i
=
0
;
i
<
tot
;
i
++
)
rev
[
i
]
=
(
rev
[
i
>>
1
]
>>
1
)
|
((
i
&
1
)
<<
(
bit
-
1
));
fft
(
a
,
1
),
fft
(
b
,
1
);
//DFT
for
(
int
i
=
0
;
i
<
tot
;
i
++
)
a
[
i
]
=
a
[
i
]
*
b
[
i
]
fft
(
a
,
-
1
);
//IDFT
for
(
int
i
=
0
;
i
<=
n
+
m
;
i
++
)
printf
(
“%d ”
,
(
int
)(
a
[
i
]。
x
/
tot
+
0。5
));
return
0
;
}
NTT
NTT程式碼(51nod1028):
#include
#include
#include
#include
#include
using
namespace
std
;
typedef
long
long
LL
;
const
LL
mod
=
998244353
;
const
int
N
=
1
<<
20
;
LL
Pow
(
LL
x
,
LL
y
){
if
(
!
y
)
return
1LL
;
LL
xx
=
Pow
(
x
,
y
/
2
);
xx
=
xx
*
xx
%
mod
;
if
(
y
&
1LL
)
xx
=
xx
*
x
%
mod
;
return
xx
;
}
LL
A
,
B
,
a
[
N
],
b
[
N
],
R
[
N
],
g
[
N
],
n
,
L
;
char
str
[
N
];
void
read
(){
scanf
(
“%s”
,
str
);
A
=
strlen
(
str
);
for
(
int
i
=
0
;
i
<
A
;
i
++
)
a
[
A
-
i
-
1
]
=
str
[
i
]
-
‘0’
;
scanf
(
“%s”
,
str
);
B
=
strlen
(
str
);
for
(
int
i
=
0
;
i
<
B
;
i
++
)
b
[
B
-
i
-
1
]
=
str
[
i
]
-
‘0’
;
}
void
NTT
(
LL
a
[
N
],
int
n
){
for
(
int
i
=
0
;
i
<
n
;
i
++
)
if
(
i
<
R
[
i
])
swap
(
a
[
i
],
a
[
R
[
i
]]);
for
(
int
d
=
1
,
t
=
n
>>
1
;
d
<
n
;
d
<<=
1
,
t
>>=
1
)
for
(
int
i
=
0
;
i
<
n
;
i
+=
(
d
<<
1
))
for
(
int
j
=
0
;
j
<
d
;
j
++
){
LL
tmp
=
g
[
t
*
j
]
*
a
[
i
+
j
+
d
]
%
mod
;
a
[
i
+
j
+
d
]
=
(
a
[
i
+
j
]
-
tmp
+
mod
)
%
mod
;
a
[
i
+
j
]
=
(
a
[
i
+
j
]
+
tmp
)
%
mod
;
}
}
int
main
(){
read
();
for
(
n
=
1
,
L
=
0
;
n
<=
A
+
B
;
n
<<=
1
,
L
++
);
for
(
int
i
=
0
;
i
<
n
;
i
++
)
R
[
i
]
=
(
R
[
i
>>
1
]
>>
1
)
|
((
i
&
1
)
<<
(
L
-
1
));
g
[
0
]
=
1
,
g
[
1
]
=
Pow
(
3
,(
mod
-
1
)
/
n
);
for
(
int
i
=
2
;
i
<
n
;
i
++
)
g
[
i
]
=
g
[
i
-
1
]
*
g
[
1
]
%
mod
;
NTT
(
a
,
n
),
NTT
(
b
,
n
);
for
(
int
i
=
0
;
i
<
n
;
i
++
)
a
[
i
]
=
a
[
i
]
*
b
[
i
]
%
mod
;
g
[
0
]
=
1
,
g
[
1
]
=
Pow
(
g
[
1
],
mod
-
2
);
for
(
int
i
=
2
;
i
<
n
;
i
++
)
g
[
i
]
=
g
[
i
-
1
]
*
g
[
1
]
%
mod
;
NTT
(
a
,
n
);
LL
Inv
=
Pow
(
n
,
mod
-
2
);
for
(
int
i
=
0
;
i
<
n
;
i
++
)
a
[
i
]
=
a
[
i
]
*
Inv
%
mod
;
for
(
int
i
=
0
;
i
<
n
-
1
;
i
++
)
a
[
i
+
1
]
+=
a
[
i
]
/
10
,
a
[
i
]
%=
10
;
int
d
;
for
(
d
=
n
-
1
;
d
&&!
a
[
d
];
d
——
);
for
(
int
i
=
d
;
i
>=
0
;
i
——
)
putchar
(
a
[
i
]
+
‘0’
);
return
0
;
}
寫在後面:博主本人時間有限(懶),對於一些基本的概念有些直接從其他部落格截圖,如果原作者看到可聯絡我刪除。