Peter_Matthew的博客

多项式

2019-02-11

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
#include<bits/stdc++.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
using namespace std;
int n,m;
struct fushu
{
double x,y;
}f[4000005],g[4000005];
fushu operator+(fushu a,fushu b){return (fushu){a.x+b.x,a.y+b.y};}
fushu operator-(fushu a,fushu b){return (fushu){a.x-b.x,a.y-b.y};}
fushu operator*(fushu a,fushu b){return (fushu){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void fft(int lim,fushu *a,int tp)
{
if(lim==1)return ;
fushu a1[lim>>1],a2[lim>>1];
for(int i=0;i<=lim;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fft(lim>>1,a1,tp);
fft(lim>>1,a2,tp);
fushu Wn=(fushu){cos(2.0*M_PI/lim),tp*sin(2.0*M_PI/lim)},w=(fushu){1.0};
for(int i=0;i<(lim>>1);i++,w=w*Wn)
{
a[i]=a1[i]+w*a2[i];
a[i+(lim>>1)]=a1[i]-w*a2[i];
}
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)
scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++)
scanf("%lf",&g[i].x);
int lim=1;
while(lim<=n+m)
lim<<=1;
fft(lim,f,1);
fft(lim,g,1);
for(int i=0;i<=lim;i++)
f[i]=f[i]*g[i];
fft(lim,f,-1);
for(int i=0;i<=n+m;i++)
printf("%.0lf ",floor(f[i].x/lim));
}

其中可以小优化一下

1
2
3
4
5
6
for(int i=0;i<(lim>>1);i++,w=w*Wn)
{
complex t=w*a2[i];
a[i]=a1[i]+t,
a[i+(lim>>1)]=a1[i]-t;
}

迭代版

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
#include<bits/stdc++.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
using namespace std;
int n,m;
struct fushu
{
double x,y;
}f[10000005],g[10000005];
fushu operator+(fushu a,fushu b){return (fushu){a.x+b.x,a.y+b.y};}
fushu operator-(fushu a,fushu b){return (fushu){a.x-b.x,a.y-b.y};}
fushu operator*(fushu a,fushu b){return (fushu){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
int l,r[10000005];
void fft(int lim,fushu *a,int tp)
{
for(int i=0;i<lim;i++)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
fushu Wn=(fushu){cos(M_PI/mid),tp*sin(M_PI/mid)};
for(int R=mid<<1,j=0;j<lim;j+=R)
{
fushu w=(fushu){1,0};
for(int k=0;k<mid;k++,w=w*Wn)
{
fushu x=a[j+k],y=w*a[j+mid+k];
a[j+k]=x+y;
a[j+mid+k]=x-y;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)
scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++)
scanf("%lf",&g[i].x);
int lim=1;
while(lim<=n+m)lim<<=1,l++;
for(int i=0;i<lim;i++)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fft(lim,f,1);
fft(lim,g,1);
for(int i=0;i<=lim;i++)
f[i]=f[i]*g[i];
fft(lim,f,-1);
for(int i=0;i<=n+m;i++)
printf("%.0lf ",floor(f[i].x/lim+0.5));
}

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
#include<bits/stdc++.h>
using namespace std;
#define int long long
const long long P=998244353,G=3,Gi=332748118;
int n,m,l,r[3000005];
long long f[3000005],g[3000005];
long long pow(long long a,long long k)
{
long long ans=1;
while(k)
{
if(k&1)ans=(ans*a)%P;
a=(a*a)%P;
k>>=1;
}
return ans%P;
}
void ntt(int lim,long long *a,int tp)
{
for(int i=0;i<lim;i++)
if(i<r[i])
swap(a[i],a[r[i]]);
for(int mid=1;mid<lim;mid<<=1)
{
long long Wn=pow(tp==1?G:Gi,(P-1)/(mid<<1));
for(int R=mid<<1,j=0;j<lim;j+=R)
{
long long w=1;
for(int k=0;k<mid;k++,w=(w*Wn)%P)
{
int x=a[j+k],y=w*a[j+mid+k]%P;
a[j+k]=(x+y)%P;
a[j+mid+k]=(x-y+P)%P;
}
}
}
}
signed main()
{
scanf("%lld%lld",&n,&m);
for(int i=0;i<=n;i++)
scanf("%lld",&f[i]),f[i]=(f[i]+P)%P;
for(int i=0;i<=m;i++)
scanf("%lld",&g[i]),g[i]=(g[i]+P)%P;
int lim=1;
while(lim<=n+m)lim<<=1,l++;
for(int i=0;i<lim;i++)
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
ntt(lim,f,1);
ntt(lim,g,1);
for(int i=0;i<=lim;i++)
f[i]=(f[i]*g[i])%P;
ntt(lim,f,-1);
long long inv=pow(lim,P-2);
for(int i=0;i<=n+m;i++)
printf("%lld ",(f[i]*inv)%P);
}
标签: 述学
使用支付宝打赏
使用微信打赏

若你觉得我的文章对你有帮助,欢迎点击上方按钮对我打赏

扫描二维码,分享此文章