概念

快速傅里叶变换(FFT)
可以使$O(n^2)$的多项式乘法运算复杂度降至$O(nlogn)$

基本原理:

  • 将多项式的系数向量转化为点值表示(DFT) $O(nlogn)$
  • 点值表示乘法 $O(n)$
  • 将多项式的点值表示转化为系数向量(IDFT) $O(nlogn)$

应用

求卷积

  • 现有两个函数$f(n)$,$g(n)$,定义$f$和$g$的卷积为$f \otimes g$
    $(f \otimes g)(n) = \sum_{i=0}^n f(i)g(n-i)$

  • 对于某些形如$h(k) = \sum_{i=0}^n f(i)g(i + k)$的问题
    可以令$f’(x) = f(n - x)$,转化为求$\sum_{i=0}^n f’(n-i)g(i + k)$

模板

DFT

复(实)数FFT,运算较快,但会有精度问题

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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>

using namespace std;

const int N = 500005;
const double pi = acos(-1.0);

char s1[N],s2[N];
int len,res[N];

struct Complex
{
double r,i;
Complex(double r=0,double i=0):r(r),i(i) {};
Complex operator+(const Complex &rhs)
{
return Complex(r + rhs.r,i + rhs.i);
}
Complex operator-(const Complex &rhs)
{
return Complex(r - rhs.r,i - rhs.i);
}
Complex operator*(const Complex &rhs)
{
return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i);
}
} va[N],vb[N];

void rader(Complex F[],int len) //len = 2^M,reverse F[i] with F[j] j为i二进制反转
{
int j = len >> 1;
for(int i = 1;i < len - 1;++i)
{
if(i < j) swap(F[i],F[j]); // reverse
int k = len>>1;
while(j>=k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}

void FFT(Complex F[],int len,int t)
{
rader(F,len);
for(int h=2;h<=len;h<<=1)
{
Complex wn(cos(t*2*pi/h),sin(t*2*pi/h));
for(int j=0;j<len;j+=h)
{
Complex E(1,0); //旋转因子
for(int k=j;k<j+h/2;++k)
{
Complex u = F[k];
Complex v = E*F[k+h/2];
F[k] = u+v;
F[k+h/2] = u-v;
E=E*wn;
}
}
}
if(t==-1) //IDFT 逆DFT
for(int i=0;i<len;++i)
F[i].r /= len;
}

void Conv(Complex a[],Complex b[],int len) //求卷积(多项式乘法)
{
FFT(a,len,1);
FFT(b,len,1);
for(int i=0;i<len;++i) a[i] = a[i]*b[i];
FFT(a,len,-1);
}

/* 上面为模板部分 */
void init(char *s1,char *s2)
{
int n1 = strlen(s1),n2 = strlen(s2);
len = 1;
while(len < 2*n1 || len < 2*n2) len <<= 1;
int i;
for(i=0;i<n1;++i)
{
va[i].r = s1[n1-i-1]-'0';
va[i].i = 0;
}
while(i<len)
{
va[i].r = va[i].i = 0;
++i;
}
for(i=0;i<n2;++i)
{
vb[i].r = s2[n2-i-1]-'0';
vb[i].i = 0;
}
while(i<len)
{
vb[i].r = vb[i].i = 0;
++i;
}
}

void gao()
{
Conv(va,vb,len);
memset(res,0,sizeof res);
for(int i=0;i<len;++i)
{
res[i]=va[i].r + 0.5;
}
for(int i=0;i<len;++i)
{
res[i+1]+=res[i]/10;
res[i]%=10;
}
int high = 0;
for(int i=len-1;i>=0;--i)
{
if(res[i])
{
high = i;
break;
}
}
for(int i=high;i>=0;--i) putchar('0'+res[i]);
puts("");
}


int main()
{
while(scanf("%s %s",s1,s2)==2)
{
init(s1,s2);
gao();
}
return 0;
}

NTT

快速数论变换,要求模数为$a2^k + 1$形式的素数

  • NTT求大整数乘法
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include<cstdio>
#include<cstring>
#include<algorithm>

using namespace std;

const int N = 400005;
const int G = 3;
const int NUM = 21;
const long long mod = (479<<21)+1;

typedef long long LL;
LL wn[30];
int len;
char s1[100010],s2[100010];
LL a[N],b[N];

LL q_pow(LL x,LL y,LL p)
{
LL res = 1;
while(y>0)
{
if(y&1) res=res*x%p;
x=x*x%p;
y>>=1;
}
return res;
}

void Getwn()
{
for(int i=0;i<21;++i)
{
int t = 1<<i;
wn[i] = q_pow(G,(mod-1)/t,mod);
}
}

void rader(LL F[],int len)
{
int j = len>>1;
for(int i=1;i<len-1;++i)
{
if(i < j) swap(F[i],F[j]);
int k=len>>1;
while(j>=k)
{
j -= k;
k >>= 1;
}
if(j<k) j += k;
}
}

void NTT(LL F[],int len,int t)
{
int id = 0;
rader(F,len);
for(int h=2;h<=len;h<<=1)
{
++id;
for(int j=0;j<len;j+=h)
{
LL E = 1;
for(int k = j;k<j+h/2;++k)
{
LL u = F[k];
LL v = (E*F[k+h/2])%mod;
F[k] = (u+v)%mod;
F[k+h/2] = ((u-v)%mod+mod)%mod;
E=(E*wn[id])%mod;
}
}
}

if(t == -1)
{
for(int i=1;i<len/2;++i) swap(F[i],F[len-i]);
LL inv = q_pow(len,mod-2,mod); //逆元
for(int i=0;i<len;++i) F[i] = (F[i]%mod * inv)%mod;
}
}

void Conv(LL a[],LL b[],int len)
{
NTT(a,len,1);
NTT(b,len,1);
for(int i=0;i<len;++i)
a[i] = (a[i]*b[i]) % mod;
NTT(a,len,-1);
}

void init(char *s1,char *s2)
{
int n1 = strlen(s1),n2 = strlen(s2);
len = 1;
while(len < 2*n1 || len < 2*n2) len <<= 1;
int i;
for(i=0;i<n1;++i)
{
a[i] = s1[n1-i-1]-'0';
}
while(i<len)
{
a[i] = 0;
++i;
}
for(i=0;i<n2;++i)
{
b[i] = s2[n2-i-1]-'0';
}
while(i<len)
{
b[i] = 0;
++i;
}
}

void gao()
{
Conv(a,b,len);
for(int i=0;i<len;++i)
{
a[i+1]+=a[i]/10;
a[i]%=10;
}
int high = 0;
for(int i=len-1;i>=0;--i)
{
if(a[i])
{
high = i;
break;
}
}
for(int i=high;i>=0;--i) putchar('0'+a[i]);
puts("");
}


int main()
{
Getwn(); //注意初始化
while(scanf("%s %s",s1,s2)==2)
{
init(s1,s2);
gao();
}
return 0;
}

题目

HDU1402

大整数乘法

HDU5307

令si表示前缀和,构造一个多项式:
$$(\sum ix^{s_i})(\sum x^{-s_i-1}) - (\sum x^{s_i})(\sum (i-1)x^{-s_i-1})$$
注意第0项要另外处理
用NTT会超时(不过标程姿势比较优雅,可以过),用DFT有精度问题,不过开long double可以水过

HDU4609

求从一些数字中选出三个能构成三角形的概率,方法是FFT求出两个构成的方案然后处理
共有$C_n^3$种取法,只需求出能构成三角形的情况,假设三角形三边长为a,b,c且a最长,则每个三角形满足

$ a > b + c , a > b , a > c$

我们可以由FFT求出所有的 b+c
然后处理成前缀和,满足第一个条件的个数即为 cnt = sum[len] - sum[a] (len 为总长度)
对每个a进行处理,将长度排序,假设a为第i个
然而有不满足其余两个条件的,现在开始去掉一些

  • $a \le b , a \le c$ : cnt -= (n-i-1)*(n-i-2)/2
  • $a \le b , a > c$ : cnt -= (n-i-1)*i
  • b + c中一条为a : cnt -= n-1

现在就找到了所有符合条件的情况

参考

  • FFT(基本上是《算法导论》的内容)
  • NTT