# FFT,NTT

20180303

https://www.cnblogs.com/rvalue/p/7351400.html (资料1)

http://blog.csdn.net/ACdreamers/article/details/39005227

http://blog.csdn.net/leo_h1104/article/details/51615710

http://blog.csdn.net/tt2767/article/details/47301849

http://blog.csdn.net/Tag_king/article/details/46351821

20190118

（复数作为$(xy)^z = x^zy^z$中的x,y, $(x^y)^z = x^{yz}$中的x, $x^yx^z = x^{y+z}$中的x时，这些式子不一定成立）

（仅非0复数的整数幂函数满足指数运算律，当指数为分数时结果是”多值“的，参考http://blog.sciencenet.cn/blog-826653-900633.html

n次单位复根指满足$\omega ^n=1$的复数$\omega$。

1.引理(消去引理)

2.引理(??)

$\omega _n^m=-\omega _n^{m+\frac{n}{2}}$

3.引理(求和引理)

（以下均假设n已经被补成2的幂）

DFT

IDFT

FFT实现

 1 #include<cstdio>
2 #include<algorithm>
3 #include<cstring>
4 #include<vector>
5 #include<cmath>
6 using namespace std;
7 #define fi first
8 #define se second
9 #define mp make_pair
10 #define pb push_back
11 typedef long long ll;
12 typedef unsigned long long ull;
13 typedef pair<int,int> pii;
14 struct cpl
15 {
16     double x,y;
17     cpl(double x=0,double y=0):x(x),y(y){}
18 };
19 cpl operator+(const cpl &a,const cpl &b)
20 {
21     return (cpl){a.x+b.x,a.y+b.y};
22 }
23 cpl operator-(const cpl &a,const cpl &b)
24 {
25     return (cpl){a.x-b.x,a.y-b.y};
26 }
27 cpl operator*(const cpl &a,const cpl &b)
28 {
29     return (cpl){a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};
30 }
31 const double pi=acos(-1);
32 const int N=2097152;
33 int rev[N];
34 void init(int len)
35 {
36     int bit=0,i;
37     while((1<<(bit+1))<=len)    ++bit;
38     for(i=0;i<len;++i)
39         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
40 }
41 void dft(cpl *a,int len,int idx)//要求len为2的幂
42 {
43     int i,j,k;cpl t1,t2,wn,wnk;
44     for(i=0;i<len;++i)
45         if(i<rev[i])
46             swap(a[i],a[rev[i]]);
47     for(i=1;i<len;i<<=1)
48     {
49         wn=cpl(cos(pi/i),idx*sin(pi/i));
50         //合并[j,j+i)和[j+i,j+2i),注意wn的n是2i而不是i
51         for(j=0;j<len;j+=(i<<1))
52         {
53             wnk=cpl(1,0);
54             for(k=j;k<j+i;++k,wnk=wnk*wn)
55             {
56                 t1=a[k];t2=a[k+i]*wnk;
57                 a[k]=t1+t2;a[k+i]=t1-t2;
58             }
59         }
60     }
61     if(idx==-1)
62     {
63         for(i=0;i<len;++i)
64             a[i].x/=len,a[i].y/=len;
65     }
66 }
67 cpl a[N],b[N];
68 int n,m;
69 int main()
70 {
71     int i;
72     init(N);
73     scanf("%d%d",&n,&m);
74     for(i=0;i<=n;++i)
75         scanf("%lf",&a[i].x);
76     for(i=0;i<=m;++i)
77         scanf("%lf",&b[i].x);
78     dft(a,N,1);dft(b,N,1);
79     for(i=0;i<N;++i)
80         a[i]=a[i]*b[i];
81     dft(a,N,-1);
82     for(i=0;i<=n+m;++i)
83         printf("%d ",int(a[i].x+0.5));
84     return 0;
85 }
View Code

NTT

（以下默认n是2的幂）（NTT部分的运算全部默认为模p意义下运算）

https://math.stackexchange.com/questions/353741/how-to-derive-this-expression-r-p-1-2-equiv-1-pmod-p-for-primitive-r/353747

https://blog.csdn.net/feynman1999/article/details/82117243

 $r*2^k+1$ r k 原根g 469762049 7 26 3 998244353 119 23 3 1004535809 479 21 3

 1 #prag\
2 ma GCC optimize(2)
3 #include<cstdio>
4 #include<algorithm>
5 #include<cstring>
6 #include<vector>
7 #include<cmath>
8 using namespace std;
9 #define fi first
10 #define se second
11 #define mp make_pair
12 #define pb push_back
13 typedef long long ll;
14 typedef unsigned long long ull;
15 const int md=998244353;
16 const int N=2097152;
17 int rev[N];
18 void init(int len)
19 {
20     int bit=0,i;
21     while((1<<(bit+1))<=len)    ++bit;
22     for(i=0;i<len;++i)
23         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
24 }
25 ll poww(ll a,ll b)
26 {
27     ll base=a,ans=1;
28     for(;b;b>>=1,base=base*base%md)
29         if(b&1)
30             ans=ans*base%md;
31     return ans;
32 }
33 void dft(int *a,int len,int idx)//要求len为2的幂
34 {
35     int i,j,k,t1,t2;ll wn,wnk;
36     for(i=0;i<len;++i)
37         if(i<rev[i])
38             swap(a[i],a[rev[i]]);
39     for(i=1;i<len;i<<=1)
40     {
41         wn=poww(idx==1?3:332748118,(md-1)/(i<<1));
42         for(j=0;j<len;j+=(i<<1))
43         {
44             wnk=1;
45             for(k=j;k<j+i;++k,wnk=wnk*wn%md)
46             {
47                 t1=a[k];t2=a[k+i]*wnk%md;
48                 a[k]+=t2;
49                 (a[k]>=md) && (a[k]-=md);
50                 a[k+i]=t1-t2;
51                 (a[k+i]<0) && (a[k+i]+=md);
52             }
53         }
54     }
55     if(idx==-1)
56     {
57         ll ilen=poww(len,md-2);
58         for(i=0;i<len;++i)
59             a[i]=a[i]*ilen%md;
60     }
61 }
62 int a[N],b[N];
63 int n,m;
64 int main()
65 {
66     int i;
67     init(N);
68     scanf("%d%d",&n,&m);
69     for(i=0;i<=n;++i)
70         scanf("%d",&a[i]);
71     for(i=0;i<=m;++i)
72         scanf("%d",&b[i]);
73     dft(a,N,1);dft(b,N,1);
74     for(i=0;i<N;++i)
75         a[i]=ll(a[i])*b[i]%md;
76     dft(a,N,-1);
77     for(i=0;i<=n+m;++i)
78         printf("%d ",a[i]);
79     return 0;
80 }
View Code

 1 #include<cstdio>
2 #include<algorithm>
3 #include<cstring>
4 #include<vector>
5 #include<cmath>
6 using namespace std;
7 #define fi first
8 #define se second
9 #define mp make_pair
10 #define pb push_back
11 typedef long long ll;
12 typedef unsigned long long ull;
13 const ll md=998244353;
14 const int bit=21,N=2097152,invN=998243877;
15 int rev[N];
16 ll wnk[21][N],iwnk[21][N];
17 //wnk[i][j]:w_{2^(i+1)}^j;iwnk[i][j]:inv(wnk[i][j])
18 //w_{2^(i+1)}=3^{(md-1)/(2^(i+1))}
19 ll poww(ll a,ll b)
20 {
21     ll base=a,ans=1;
22     for(;b;b>>=1,base=base*base%md)
23         if(b&1)
24             ans=ans*base%md;
25     return ans;
26 }
27 void init(int len)
28 {
29     int i,j,ed;ll wn;
30     for(i=0;i<len;++i)
31         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
32     for(i=0;i<bit;++i)
33     {
34         ed=1<<(i+1);
35         wnk[i][0]=iwnk[i][0]=1;
36         wn=wnk[i][1]=poww(3,(md-1)/(1<<(i+1)));
37         for(j=2;j<ed;++j)
38             wnk[i][j]=wnk[i][j-1]*wn%md;
39         wn=iwnk[i][1]=poww(332748118,(md-1)/(1<<(i+1)));
40         for(j=2;j<ed;++j)
41             iwnk[i][j]=iwnk[i][j-1]*wn%md;
42     }
43 }
44 void dft(ll *a,int len,int idx)//要求len为2的幂
45 {
46     int i,i1,j,k;ll t1,t2;
47     for(i=0;i<len;++i)
48         if(i<rev[i])
49             swap(a[i],a[rev[i]]);
50     for(i=1,i1=0;i<len;i<<=1,++i1)
51     {
52         const ll *wnk=idx==1?(::wnk[i1]):iwnk[i1];
53         for(j=0;j<len;j+=(i<<1))
54         {
55             for(k=j;k<j+i;++k)
56             {
57                 t1=a[k];t2=a[k+i]*wnk[k-j]%md;
58                 a[k]+=t2;a[k+i]=t1-t2;
59                 (a[k]>=md) && (a[k]-=md);
60                 (a[k+i]<0) && (a[k+i]+=md);
61             }
62         }
63     }
64     if(idx==-1)
65     {
66         for(i=0;i<len;++i)
67             (a[i]*=invN)%=md;
68     }
69 }
70 ll a[N],b[N];
71 int n,m;
72 int main()
73 {
74     int i;
75     init(N);
76     scanf("%d%d",&n,&m);
77     for(i=0;i<=n;++i)
78         scanf("%lld",&a[i]);
79     for(i=0;i<=m;++i)
80         scanf("%lld",&b[i]);
81     dft(a,N,1);dft(b,N,1);
82     for(i=0;i<N;++i)
83         a[i]=a[i]*b[i]%md;
84     dft(a,N,-1);
85     for(i=0;i<=n+m;++i)
86         printf("%lld ",a[i]);
87     return 0;
88 }
View Code

todo

posted @ 2019-01-27 22:19  hehe_54321  阅读(284)  评论(0编辑  收藏  举报