FFT小结

Posted 2020-10-25 LadyLex

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了FFT小结相关的知识，希望对你有一定的参考价值。

FFT小结！

零，说在前面

……令人毒瘤的FFT

其实FFT是卷积这个庞大学术课题的一个小应用啦.....

WC的时候听卷积定理真是上天的体验

那么话不多说我们开始FFT的总结.

板子我们就不放了!这个是很基础的,网上也有讲解的很详细的学习资料

大家请善用搜索引擎啦......

一，FFT与NTT的简单例题以及基础思想&技巧

首先要提到的就是反转思想了。

可能有的题目形式并不是卷积，但是我们通过改变枚举顺序，或者反转数组就可以得到卷积的形式

这样的一个典型题目就是bzoj2194了，非常基础的题目……

这里就不推式子了

然后是关于字符串匹配的问题……说一种最简单的做法吧，我们把一个串反过来，然后同种字符设为1，其余的为0，然后做卷积

这样就可以拿到哪些位置匹配了，然后把每种字符都做一遍就可以得出2个字符串有多长是一样的……

的确有的题问这个的……

以及FFT还可以和斯特林数结合起来……然后这个我正在学+做题，做完再来补上

二，FFT的更多应用

其实我和Cooook写了一个课件的包括了底下大部分东西……啊我好懒啊我不想再打一遍了

如果有想要课件的同学跟我说下啊……我把密码放在评论区……

链接

一，多项式求逆

二，多项式开根

三，多项式除法&取模

四，多项式取$ln$

五，多项式取$exp$

六，多项式快速插值

七，多项式多点求值

八，指数型母函数

九，任意模数FFT

这个我还没有看……现在正在学ing，请等我补坑……

十，与卷积定理的结合

这个在WC2018上听LCA爷爷讲了讲……听傻了……

现在我还是不会的…估计一段时间内我也不太可能会…留坑待填

三，FFT/NTT好题选讲

后面的题我在不断做不断更新中……

一，UOJ#50

这个题我们发现其实就是一个类似于二叉树计数的题目。

首先我们发现那个标号的限制完全没有用……

不难想到一个$O(n^{3})$的dp方程，即枚举两个中子打到的子树大小，然后再看剩下的大小是否符合破坏死光的要求，再乘组合数。这是一个40pts做法。

然后，这个dp出来之后，我们发现儿子的贡献就是一个卷积的形式啊！

然后把合并写成指数型母函数的卷积，这样的复杂度是$O(n^{2}log_{n})$,可以拿到60pts

然后……你会很快的发现，如果我们把死光也看成多项式，那么就是3个卷起来，那么复杂度就成了$O(nlog_{n}^{2})$，用CDQ+FFT可以得到90pts～100pts

那么我个人就是打的这种方法！4s打2e5的CDQ+FFT是没有问题的！

最有意思的就是那个多项式的维护了。由于我们的式子是这样的：

$$F‘(x)=C(x)*F^{2}(x)+1$$

所以我们要维护一个平方的多项式，代表表示1～区间中点mid平方的结果

每次我们用平方数组更新f，然后用新计算了mid+1~r，我们就把它更新到原来那个平方里面去

怎么更新？提示：$(x+a)^{2}=x^{2}+2ax+a^{2}$

当然，我们还有一种解微分方程的做法，但是鉴于我的数学水平我并不会23333。

有兴趣的同学可以看Vfk的题解。

给个代码：

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <cstdlib>
 5 using namespace std;
 6 #define N 400010
 7 #define L 1100000
 8 #define mod 998244353
 9 #define phi 998244352
10 #define g 3
11 #define LL long long
12 inline int max(int a,int b){return a>b?a:b;}
13 inline int min(int a,int b){return a<b?a:b;}
14 char s[N];
15 int f[N],pf[N],n,logi[N],litlen=1,len=1;
16 int bin[25],poww[L],rev[L],fac[L],inv[L],invf[L];
17 inline int quick_mod(int di,int mi)
18 {
19     int ret=1;
20     for(;mi;mi>>=1,di=(LL)di*di%mod)
21         if(mi&1)ret=(LL)ret*di%mod;
22     return ret;
23 }
24 inline void upd(int &a,int b){a+=b;if(a>=mod)a-=mod;}
25 inline int sub(int a,int b){a-=b;if(a<0)a+=mod;return a;}
26 inline int dft(int *a,int opt,int le)
27 {
28     register int i,j,d,l,wn,w,tmp;
29     for(i=0;i<le;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
30     for(d=2;d<=le;d<<=1)
31         for(wn=(opt==1)?poww[len/d]:poww[len-len/d],l=d>>1,i=0;i<le;i+=d)
32             for(w=1,j=0;j<l;++j,w=(LL)w*wn%mod)
33                 tmp=(LL)a[i+j+l]*w%mod,a[i+j+l]=(a[i+j]-tmp+mod)%mod,a[i+j]=(a[i+j]+tmp)%mod;
34                 // tmp=(LL)a[i+j+l]*w%mod,a[i+j+l]=sub(a[i+j],tmp),upd(a[i+j],tmp);
35     if(opt==-1)
36         for(i=0;i<le;++i)a[i]=(LL)a[i]*inv[le]%mod;
37 }
38 int ret[L],C[L],A[L],B[L];
39 inline void CDQ(int l,int r)
40 {
41     if(l==r)
42         if(l<n){printf("%d\n",(LL)f[l]*fac[l-1]%mod);return;}
43         else {printf("%d\n",(LL)f[l]*fac[l-1]%mod);exit (0);}
44     register int le=r-l+1,i,mi=(l+r)>>1;
45     CDQ(l,mi);
46     for(i=0;i<le;++i)
47         if(i&1)rev[i]=(rev[i>>1]>>1)|(le>>1);
48         else rev[i]=(rev[i>>1])>>1;
49     memset(A,0,le<<2),memset(B,0,le<<2);
50     for(i=l;i<=mi;++i)A[i-l]=pf[i];
51     for(i=0;i<r-l;++i)B[i]=C[i];
52     dft(A,1,le),dft(B,1,le);
53     for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod;
54     dft(ret,-1,le);
55     for(i=mi+1;i<=r;++i)upd(f[i],ret[i-l-1]);
56     memset(A,0,le<<2),memset(B,0,le<<2);
57     for(i=l;i<=mi;++i)A[i-l]=(LL)f[i]*inv[i]%mod;
58     for(i=0;i<=r-l;++i)
59         if(i<l)B[i]=(LL)f[i]*inv[i]*2%mod;
60         else if(i<=mi)B[i]=(LL)f[i]*inv[i]%mod;
61     dft(A,1,le),dft(B,1,le);
62     for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod;
63     dft(ret,-1,le);
64     for(i=mi+1;i<=r;++i)upd(pf[i],ret[i-l]);
65     CDQ(mi+1,r);
66 }
67 int main()
68 {
69     register int i,j;
70     for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1;
71     scanf("%d%s",&n,s);
72     while(len<(n<<1))len<<=1;
73     poww[0]=1,poww[1]=quick_mod(g,phi/len);
74     for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
75     for(fac[0]=fac[1]=1,i=2;i<=len;++i)fac[i]=(LL)fac[i-1]*i%mod;
76     for(inv[0]=inv[1]=1,i=2;i<=len;++i)inv[i]=(LL)inv[mod%i]*(mod-mod/i)%mod;
77     for(invf[0]=invf[1]=1,i=2;i<=len;++i)invf[i]=(LL)invf[i-1]*inv[i]%mod;
78     for(i=0;i<n;++i)
79         if(s[i]==‘1‘)C[i]=(LL)invf[i]*inv[2]%mod;
80     f[1]=1,CDQ(1,len>>1);
81 }

UOJ#50

二，UOJ#182

这题我似乎在课件里面讲了！

然后写的自己还算满意吧……再简单说一下

总之你可以发现每个数再若干次操作之后可以写成$\frac {Ax+C}{Bx+D}$的形式，x代表每个数的初始值。然后……你把它写丑一点，可以变成形如$e+f*\frac{1}{x+g}$的形式……但是具体内容请自己推导吧……

那么总和就是$n*e+f*\sum\frac{1}{x+g}$，问题转化为维护$\sum\frac{1}{x_{i}+g}$

如果我们把g看成变量的话

那么这个东西，通分可以得到等价于$\frac{ \sum _{j} \prod _{i!=j} {A_{i}+x} } {\prod _{j} (A_{j}+x)}$

注意到分母里的式子（设为$f$）可以倍增加FFT求出来

分子的式子（设为$g$），经过我们的推导，是$f$的导数

你问我怎么推导？

第一种方法，你可以打表

第二种方法，观察特殊点……我当时只看到了$[x^{n-1}]g(x)=n$

第三种方法，我们推他一下……

啊……我推了半小时推不出来

如果有会的同学请在评论区教教我

然后呢？我们可以套用多项式多点求值的模板，求出两个函数每一个点的点值，然后计算答案就行了

然后这题让我知道了指针的常数很小23333

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <ctime>
  4 #include <algorithm>
  5 using namespace std;
  6 #define mod 998244353
  7 #define LL long long
  8 #define RG register
  9 #define N 100010
 10 #define M 60010
 11 char cB[1<<15],*S=cB,*T=cB;
 12 #define getc (S==T&&(T=(S=cB)+fread(cB,1,1<<15,stdin),S==T)?0:*S++)
 13 inline int read()
 14 {
 15     int x=0;RG char c=getc;
 16     while(c<‘0‘|c>‘9‘)c=getc;
 17     while(c>=‘0‘&c<=‘9‘)x=10*x+(c^48),c=getc;
 18     return x;
 19 }
 20 inline int quick_mod(int di,int mi)
 21 {
 22     int ret=1;
 23     for(;mi;mi>>=1,di=(LL)di*di%mod)
 24         if(mi&1)ret=(LL)ret*di%mod;
 25     return ret;
 26 }
 27 #define L (1<<18)+10
 28 #define LM (1<<16)+10
 29 int logi[L],inv[L],bin[25],rev[L],poww[L],len=1,ro=1;
 30 inline void dft(int *a,int n,int opt)
 31 {
 32     RG int i,j,d=logi[len]-logi[n],wn,l,tmp,*x,*y,*w;
 33     for(i=0;i<n;++i)if(i<(rev[i]>>d))swap(a[i],a[rev[i]>>d]);
 34     for(d=2;d<=n;d<<=1)
 35         for(i=0,l=d>>1,wn=(opt==1)?(len/d):(-len/d);i<n;i+=d)
 36             for(w=poww+((opt==1)?0:len),j=0,x=a+i+j,y=x+l;j<l;++j,w+=wn,++x,++y)
 37                 tmp=(LL)(*w)*(*y)%mod,*y=(*x-tmp+mod)%mod,*x=(*x+tmp)%mod;
 38     if(opt==-1)for(i=0;i<n;++i)a[i]=(LL)a[i]*inv[n]%mod;
 39 }
 40 int n,m,cnt,sum,x[N];
 41 int g[L<<1],h[L<<1];
 42 int sta[LM],ans[M];
 43 struct node{int e,f,g,id;}s[M];
 44 inline void ins(int a,int b,int c,int d,int id)
 45 {
 46     if(b==0)
 47         {ans[id]=(((LL)a*sum)+((LL)c*n))%mod*(LL)quick_mod(d,mod-2)%mod;return;}
 48     s[++cnt].f=((LL)b*c-(LL)a*d)%mod,
 49     b=quick_mod(b,mod-2);
 50     s[cnt].e=(LL)a*b%mod,s[cnt].g=(LL)d*b%mod,
 51     b=(LL)b*b%mod,s[cnt].f=(LL)s[cnt].f*b%mod;
 52     if(s[cnt].f<0)s[cnt].f+=mod;
 53     s[cnt].id=id;
 54 }
 55 inline int solve0(int *a,int l,int r)
 56 {
 57     RG int ra=r-l+2;
 58     if(ra<=128)
 59     {
 60         memset(a,0,ra<<2),a[0]=1;
 61         RG int i,j,v;
 62         for(i=l;i<=r;++i)
 63             for(v=x[i],j=i-l;~j;--j)
 64                 a[j+1]=(a[j+1]+a[j])%mod,
 65                 a[j]=(LL)a[j]*v%mod;
 66         return ra;
 67     }
 68     RG int i,mi=l+r>>1,la=1;
 69     while(la<ra)la<<=1;
 70     RG int *f1=a,r1=solve0(f1,l,mi);
 71     RG int *f2=a+r1,r2=solve0(f2,mi+1,r);
 72     static int tmp1[L],tmp2[L];
 73     memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1);
 74     memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1);
 75     for(i=0;i<la;++i)a[i]=(LL)tmp1[i]*tmp2[i]%mod;
 76     dft(a,la,-1);
 77     return ra;
 78 }
 79 int *p[L];
 80 inline int solve1(int id,int l,int r)
 81 {
 82     static int mem[LM<<4],*head=mem;
 83     RG int ra=r-l+2;
 84     if(ra<=128)
 85     {
 86         RG int i,j,v,*f=p[id]=head;head+=ra;
 87         memset(f,0,ra<<2),f[0]=1;
 88         for(i=l;i<=r;++i)
 89             for(v=(mod-sta[i])%mod,j=i-l;~j;--j)
 90                 f[j+1]=(f[j+1]+f[j])%mod,f[j]=(LL)f[j]*v%mod;
 91         return ra;
 92     }
 93     RG int mi=l+r>>1,la=1;
 94     while(la<ra)la<<=1;
 95     RG int r1=solve1(id<<1,l,mi),*f1=p[id<<1];
 96     RG int r2=solve1(id<<1|1,mi+1,r),*f2=p[id<<1|1];
 97     static int tmp1[LM],tmp2[LM];
 98     memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1);
 99     memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1);
100     RG int i,*f=p[id]=head;head+=ra;
101     for(i=0;i<la;++i)f[i]=(LL)tmp1[i]*tmp2[i]%mod;
102     dft(f,la,-1);return ra;
103 }
104 inline int get_inv(int *a,int *ret,int ra)
105 {
106     if(ra==1){ret[0]=quick_mod(a[0],mod-2);return 1;}
107     static int tmp[L];
108     RG int r1=get_inv(a,ret,(ra+1)>>1),i,la=1;
109     while(la<(ra<<1))la<<=1;
110     memcpy(tmp,a,ra<<2),memset(tmp+ra,0,(la-ra)<<2);
111     memset(ret+r1,0,(la-r1)<<2);
112     dft(tmp,la,1),dft(ret,la,1);
113     for(i=0;i<la;++i)
114     {
115         ret[i]=((LL)ret[i]*(2ll-((LL)ret[i]*tmp[i]%mod)))%mod;
116         if(ret[i]<0)ret[i]+=mod;
117     }
118     dft(ret,la,-1);return ra;
119 }
120 inline void rev_copy(int *st,int *to,int ra)
121     {for(RG int i=0;i<ra;++i)to[i]=st[ra-i-1];}
122 inline void reverse(int *st,int ra)
123     {for(RG int t,i=0,j=ra-1;i<j;++i,--j)t=st[i],st[i]=st[j],st[j]=t;}
124 inline int get_mod(int *a,int ra,int *p,int rp,int *ret)
125 {
126     while(ra&&!a[ra-1])--ra;
127     while(rp&&!p[rp-1])--rp;
128     if(ra<rp){memcpy(ret,a,ra<<2),memset(ret+ra,0,(rp-ra)<<2);return rp;}
129     static int tmp1[L],tmp2[L];
130     RG int i,j,re=ra-rp+1,la=1;
131     while(la<(re<<1))la<<=1;
132     memset(tmp1,0,la<<2),rev_copy(p,tmp1,rp),get_inv(tmp1,tmp2,re);
133     memset(tmp2+re,0,(la-re)<<2),dft(tmp2,la,1);
134     rev_copy(a,tmp1,ra),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1);
135     for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
136     dft(tmp1,la,-1);
137     la=1;while(la<ra)la<<=1;
138     reverse(tmp1,re),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1);
139     memcpy(tmp2,p,rp<<2),memset(tmp2+rp,0,(la-rp)<<2),dft(tmp2,la,1);
140     for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
141     dft(tmp1,la,-1);
142     for(i=0;i<rp;++i)ret[i]=(a[i]-tmp1[i]+mod)%mod;
143     memset(ret+rp,0,(la-rp)<<2);
144     while(rp&&!ret[rp-1])--rp;
145     return rp;
146 }
147 int val0[LM],val[LM];
148 inline void solve2(int id,int *a,int *b,int l,int r,int deg)
149 {
150     RG int ra=r-l+2;
151     if(deg<=128)
152     {
153         RG int i,j,v,u,F,G;
154         for(i=l;i<=r;val0[i]=F,val[i]=G,++i)
155             for(u=sta[i],v=1,F=G=0,j=0;j<=deg;++j,v=(LL)v*u%mod)
156                 F=(F+(LL)v*a[j])%mod,G=(G+(LL)v*b[j])%mod;
157         return;
158     }
159     RG int i,mi=l+r>>1;
160     int r1=get_mod(a,deg,p[id],ra,a+deg);a+=deg;
161     int r2=get_mod(b,deg,p[id],ra,b+deg);b+=deg;
162     ra=min(r1,r2);
163     solve2(id<<1,a,b,l,mi,ra),solve2(id<<1|1,a,b,mi+1,r,ra);
164 }
165 int main()
166 {
167     RG int i,j,opt,v,A=1,B=0,C=0,D=1;
168     n=read(),m=read();
169     while(len<(n<<1))len<<=1;
170     for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1;
171     for(i=0;i<=20;++i)
172         inv[bin[i]]=quick_mod(bin[i],mod-2),logi[bin[i]]=i;
173     for(int i=0;i<len;++i)
174         if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1);
175         else rev[i]=(rev[i>>1]>>1);
176     poww[0]=1,poww[1]=quick_mod(3,(mod-1)/len);
177     for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
178     for(sum=0,i=1;i<=n;++i)x[i]=read(),sum=(sum+x[i])%mod;
179     for(i=1;i<=m;++i)
180         if(read()==1)v=read(),A=(A+(LL)v*B)%mod,C=(C+(LL)v*D)%mod,ins(A,B,C,D,i);
181         else swap(A,B),swap(C,D),ins(A,B,C,D,i);
182     if(cnt)
183     {
184         for(i=1;i<=cnt;++i)sta[i]=s[i].g;
185         sort(sta+1,sta+cnt+1),opt=unique(sta+1,sta+cnt+1)-sta-1;
186         for(i=1;i<=cnt;++i)
187             s[i].g=lower_bound(sta+1,sta+opt+1,s[i].g)-sta;
188         solve0(g,1,n);
189         for(h[n]=0,i=1;i<=n;++i)h[i-1]=(LL)g[i]*i%mod;
190         solve1(1,1,opt),solve2(1,g,h,1,opt,n+1);
191         for(i=1;i<=cnt;++i)
192             ans[s[i].id]=( (LL)s[i].e*n+(LL)s[i].f*val[s[i].g]%mod*quick_mod(val0[s[i].g],mod-2) )%mod;
193     }
194     for(i=1;i<=m;++i)printf("%d\n",ans[i]);
195 }