FFT:
没啥好说的吧。。
证明应该都会,写的时候记住两个点就行:
1.怎么定义复数?千万别写成
complex<double> w=(1,0);
可以自己试一下这样输出什么东西……
2.枚举len,遍历前一半,用原来的$a_{i},a_{i+len/2}$值计算新的$a_{i},a_{i+len/2}$值。
我觉得我还是解释一下这么做的原因好一点。。
我们希望得到多项式$A(x)$在$w_{n}^{1},w_{n}^{2},cdots,w_{n}^{n}$处的点值表示。
而我们有$w_{n}^{k}=w_{frac{n}{2}}^{frac{k}{2}},w_{n}^{k+frac{n}{2}}=-w_{frac{n}{2}}^{frac{k}{2}}$。
那么把偶数下标的系数和奇数下标的系数分别拆出来组成两个多项式$A_{0}(x)$和$A_{1}(x)$,则有
$A(w_{n}^{k})=A_{0}(w_{n}^{2k})+w_{n}^{k}A_{1}(w_{n}^{2k})=A_{0}(w_{frac{n}{2}}^{k})+w_{n}^{k}A_{1}(w_{frac{n}{2}}^{k})$
$A(w_{n}^{k+frac{n}{2}})=A_{0}(w_{n}^{2k})-w_{n}^{k}A_{1}(w_{n}^{2k})=A_{0}(w_{frac{n}{2}}^{k})-w_{n}^{k}A_{1}(w_{frac{n}{2}}^{k})$
那么我们可以先递归处理$A_{0}(x)$和$A_{1}(x)$,然后爆枚$k$,得到$A(x)$的$n$个点值。
这东西很像线段树,倒数第$i$层的$n$是$2^{i}$,有$frac{n}{2^{i}}$个多项式需要爆枚$2^{i}$次。
每层爆枚的复杂度是$O(n)$,一共$log{n}$层,总复杂度$O(nlog{n})$。
发现位置$i$最终会被换到$i$的二进制位反过来的位置,于是可以不用递归,直接模拟线段树上爆枚的过程即可。
强烈推荐伟大的石神给出的里程碑式的证明!
代码:
#include<bits/stdc++.h> #define maxn 3000005 #define maxm 500005 #define inf 0x7fffffff #define ll long long #define cp complex<double> #define pi acos(-1) using namespace std; cp a[maxn],b[maxn],c[maxn]; int ind[maxn]; inline int read(){ int x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline void fft(cp *t,int n,int op){ for(int i=0;i<n;i++) ind[i]=(i&1)?((ind[i>>1]>>1)|(n>>1)):(ind[i>>1]>>1); for(int i=0;i<n;i++) if(ind[i]<i) swap(t[i],t[ind[i]]); for(int len=2;len<=n;len<<=1) for(int j=0;j<n;j+=len){ cp w(1,0),p(cos(op*2*pi/len),sin(op*2*pi/len)); for(int k=0;k<(len>>1);k++,w*=p){ cp t1=t[j+k+(len>>1)],t2=t[j+k]; t[j+k+(len>>1)]=t2-t1*w,t[j+k]=t2+t1*w; } } return; } int main(){ int n=read(),m=read(); for(int i=0;i<=n;i++) a[i]=read(); for(int i=0;i<=m;i++) b[i]=read(); int len=1; while(len<=n+m) len<<=1; fft(a,len,1),fft(b,len,1); for(int i=0;i<len;i++) c[i]=a[i]*b[i]; fft(c,len,-1); for(int i=0;i<=n+m;i++) cout<<(int)(c[i].real()/len+0.5)<<" "; cout<<endl; return 0; }
NTT:
原根的性质与单位根一模一样。只是$IDFT$的时候得求个逆元。
题里有模数就用$NTT$,没有模数也尽量用$NTT$,避免出现时间精度双爆炸的惨剧。
注意$w_{n}^{k}=g^{frac{mod-1}{n}cdot k}$,证明时将k分离出来考虑即可。
我还是证一下吧。。虽然这也不算是证明了。。
根据原根的性质,有$g_{i}mod{p}$的值互不相同$(0leq i<p)$。
那么实际上$w_{n}^{k}=g^{frac{k}{n}}$,都满足首尾值相等且中间值互不相同。
但是等一下,这个$frac{k}{n}$不是整数啊?
而且指数没有同余定律,只有费马小定理$a^{p-1}equiv apmod p$。
那我们只能取$n|(mod-1)$的$n$来做了,此时$g^{frac{k}{n}}=g^{frac{k(mod-1)}{n}}$。
现在我们就理解为什么$NTT$只能用形如$2^{k}+1$的模数了。
伟大的石神在刚才那篇里也更了$NTT$!
代码:
#include<bits/stdc++.h> #define maxn 3000005 #define maxm 500005 #define inf 0x7fffffff #define ll long long #define mod 998244353 #define g 3 using namespace std; ll a[maxn],b[maxn],res[maxn],ind[maxn]; inline ll read(){ ll x=0,f=1; char c=getchar(); for(;!isdigit(c);c=getchar()) if(c=='-') f=-1; for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f; } inline ll power(ll a,ll b){ll ans=1;while(b)ans=(b&1)?ans*a%mod:ans,a=a*a%mod,b>>=1;return ans;} inline void ntt(ll *t,ll n,ll op){ for(ll i=0;i<n;i++) ind[i]=(i&1)?((ind[i>>1]>>1)|(n>>1)):(ind[i>>1]>>1); for(ll i=0;i<n;i++) if(ind[i]<i) swap(t[i],t[ind[i]]); for(ll len=1;len<=n;len<<=1){ ll p=power(g,(mod-1)/len); if(op==-1) p=power(p,mod-2); for(ll i=0;i<n;i+=len) for(ll j=i,w=1,tp;j<i+(len>>1);j++,w=w*p%mod) tp=t[j+(len>>1)]*w%mod,t[j+(len>>1)]=(t[j]-tp+mod)%mod,t[j]=(t[j]+tp)%mod; } return; } int main(){ ll n=read(),m=read(); for(ll i=0;i<=n;i++) a[i]=read(); for(ll i=0;i<=m;i++) b[i]=read(); ll len=1; while(len<=n+m) len<<=1; ntt(a,len,1),ntt(b,len,1); for(ll i=0;i<len;i++) res[i]=a[i]*b[i]; ntt(res,len,-1); ll mo=power(len,mod-2); for(ll i=0;i<=n+m;i++) cout<<res[i]*mo%mod<<" "; cout<<endl; return 0; }
(upd:数组一定要开两倍,点数少了就没法表示一个多项式)