离散傅里叶变换
定义：
是将 $\omega_n^1,\omega_n^2,…,\omega_n^n$ 带入多项式所得到的一种特殊的多项式点值表示
性质：
考虑将 $A(x)$ 的每一项按下标的奇偶分成两个多项式之和：
$A(x) = (a_0 + a_2x^2 + … + a_{n - 2}x^{n - 2}) + (a_1x + a_3x^3 + … + a_{n-1}x^{n-1})$
如果设：
$A_1(x) = a_0 + a_2x + … + a_{n - 2}x^{\frac{n}{2} - 1}$
$A_2(x) = a_1 + a_3x + … + a_{n - 1}x^{\frac{n}{2} - 1}$
那么：
$A(x) = A_1(x^2) + xA_2(x^2)$
于是我们把 $x = \omega_n^k$ 带入$A(x)$：
当 $k\subseteq \lbrace 0,1,2,… ,\frac{n}{2}-1\rbrace $ 时：
$A(\omega_{n}^{k})=A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$
$A(\omega_{n}^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$
这样如果我们知道了 $A_1(x),A_2(x)$ 的离散傅里叶变换表示就可以 $O(n)$ 地求出 $A(x)$ 的离散傅里叶变换。
快速傅里叶变换
根据上面所提的离散傅里叶变换的性质，我们就可以递归地在 $O(nlogn)$的时间内将多项式的系数表示转换为点值表示。下面就介绍一下优化和实现上的细节：
蝴蝶操作：
有一个很神奇的性质：将一个多项式按下标的奇偶分为两部分，并将奇数部分放在偶数部分之后，然后不断递归下去。那么多项式下标为 $i$ 的数在这个操作之后的下标就变成了「将 $i$ 二进制翻转所得的数」。
1
| for(rint i=0; i<lim; ++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(l-1));
|
根据这个性质，我们可以先预处理出多项式中每个数操作之后的位置，并将它放到这个位置。然后类似于递归时回溯的操作，不断向上还原求出点值表示。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| for(rint i=0; i<lim; ++i)
if(i<pos[i]) swap(a[i], a[pos[i]]);
for(rint i=1; i<lim; i<<=1)
{
Complex del(cos(PI/i), sin(PI/i));
for(rint j=0; j<lim; j+=(i<<1))
{
Complex t(1, 0);
for(rint k=0; k<i; ++k, t=t*del)
{
Complex x=a[j+k], y=t*a[j+k+i];
a[j+k]=x+y; a[j+k+i]=x-y;
}
}
}
|
傅里叶逆变换
按照上述方法我们在 $O(nlogn)$ 内将多项式的系数表示转换成了点值表示，然后我们在 $O(n)$ 的时间内将两个多项式的点值表示相乘，得到了乘积的点值表示。利用傅里叶逆变换就可以同样在 $O(nlogn)$ 内将点值表示转换成系数表示。
离散傅里叶变换的另一个性质：把多项式 $A(x)$ 的离散傅里叶变换作为另一个多项式 $B(x)$ 的系数，并将 $\omega_{n}^{0}, \omega_{n}^{-1}, \omega_{n}^{-2}, …, \omega_{n}^{-(n - 1)}$ 代入多项式 $B(x)$ ，得到的每个数除 $n$，这样的 $n$ 个数就恰好是多项式 $A(x)$ 的系数。
根据上面的性质，我们在求离散傅里叶变换的代码上稍作修改，就能在 $O(nlogn)$ 内将点值表示转换成系数表示。
最终代码
下面是用 FFT 优化的高精度乘法代码（BZOJ2179）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
| #include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <vector>
#include <map>
#include <set>
#define MAXN 200005
#define INF 0x3f3f3f3f
#define rint register int
#define LL long long
#define LD long double
using namespace std;
const double PI=acos(-1);
int n, m, len, pos[MAXN], ans[MAXN];
char sa[MAXN], sb[MAXN];
struct CP
{
double x, y;
CP(double a=0, double b=0)
{
x=a; y=b;
}
friend CP operator * (CP a, CP b)
{
return CP(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
}
friend CP operator + (CP a, CP b)
{
return CP(a.x+b.x, a.y+b.y);
}
friend CP operator - (CP a, CP b)
{
return CP(a.x-b.x, a.y-b.y);
}
} a[MAXN], b[MAXN];
void fft(CP a[], int op)
{
for(rint i=0; i<n; ++i)
if(i<pos[i]) swap(a[i], a[pos[i]]);
for(rint i=1; i<n; i<<=1)
{
CP omg(cos(PI/i), op*sin(PI/i));
for(rint j=0; j<n; j+=(i<<1))
{
CP t(1, 0);
for(rint k=0; k<i; ++k, t=t*omg)
{
CP x=a[j+k], y=t*a[j+k+i];
a[j+k]=x+y; a[j+k+i]=x-y;
}
}
}
}
int main()
{
scanf("%d", &len);
scanf("%s%s", sa, sb);
for(rint i=0; i<len; ++i)
{
a[i].x=sa[len-1-i]-'0';
b[i].x=sb[len-1-i]-'0';
}
for(n=1; n<(len*2); ) n<<=1, m++;
for(rint i=0; i<=n; ++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(m-1));
fft(a, 1);
fft(b, 1);
for(rint i=0; i<n; ++i) a[i]=a[i]*b[i];
fft(a, -1);
for(rint i=0; i<n; ++i)
{
ans[i]+=(int)(a[i].x/n+0.5);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
len=len*2-1;
while(!ans[len] && len>=1) len--;
for(rint i=len; i>=0; i--) printf("%d", ans[i]);
return 0;
}
|
例题
1.BZOJ3771 Triple 2.BZOJ4259 残缺的字符串 3.AGC005F Many Easy Problems