今天在北京某培训营中dyf给我们普及了一些FFT的基本应用...不过唯一的收获是终于大致搞清楚了怎么做多项式的牛顿迭代...这样也就可以向策爷学一点新的东西了。
还是从一道题开始。
【某模拟题T2】求所有循环长度均在集合[tex]U[/tex]中的长度为[tex]n[/tex]的置换个数。[tex]n \leq 10^5[/tex]。出题人：jcvb
首先可以写出递推式：
[tex]F_n=\sum_{i \in U} F_{n-i} \times \binom{n-1}{i-1} \times (i-1)![/tex]
这显然可以用分治FFT搞，时间复杂度[tex]O(n \log^2 n)[/tex]。当然这不是重点，接下来主要讨论这道题的[tex]O(n\log n)[/tex]的解法。
令[tex]A_{(z)}[/tex]为只有一个循环的合法置换的方案数的指数生成函数，那么有：
[tex]A_{(z)}=\sum_{i \in U} \frac{1}{i} z^i[/tex]
又令[tex]F_{(z)}[/tex]为长度为最终答案的指数生成函数，那么就有：
[tex]F_{(z)}=\sum_{i=0}^{\infty} \frac{A_{(z)}^i}{i!}=e^{A_{(z)}}[/tex]
对两边求自然对数得 [tex]\ln F_{(z)}-A_{(z)}=0[/tex]
把多项式[tex]F_{(z)}[/tex]看做未知数[tex]x[/tex]，用牛顿迭代解这个方程：
[tex]x_{i+1}=x_{i}-x_{i} \times (\ln x_i-A_{(z)})[/tex]
因为牛顿迭代精度是以平方级别增加的，所以第[tex]i[/tex]次迭代后得到的前[tex]2^i[/tex]位是精确的。在迭代[tex]O(\log n)[/tex]次之后便可以得到答案。
唯一留下的问题是怎么对一个多项式求自然对数。设[tex]B_{(z)}= \ln A_{(z)}[/tex]
对等式两边求导得[tex]B_{(z)}'=\frac{A_{(z)}'}{A_{(z)}}[/tex]
而由于牛顿迭代的精度是[tex]2^i[/tex]位，所以可以再模[tex]x^{2^i}[/tex]下进行，用多项式求逆搞然后再积分回去就可以了。
最后分析复杂度我们可以得到这个算法的复杂度是[tex]O(n\log n)[/tex]的，实际上跑得的确比分治FFT快了那么一点点。
用这个算法就可以求得一个常数项为0的多项式的exp了。
【应用】多项式的快速幂（模[tex]x^n[/tex]下）
显然这个是可以用快速幂[tex]O(n \log^2 n)[/tex]搞，但是用多项式求exp就可以得到[tex]O(n \log n)[/tex]的解法。
令[tex]A_{(z)}=B_{(z)}^k[/tex]，对两边取自然对数得[tex]\ln A_{(z)}=k \ln B_{(z)}[/tex]，然后再对左边求exp就好了。
这个东西可以和一个神奇的定理形成一个强力的组合拳，但是鉴于那玩意要拿来出题，那我就不写了。
然后吐槽一下多项式算法的常数，FFT常数够大了吧？套上一层倍增求逆我还是[tex]O(n \log n)[/tex]的。多项式求逆常数够大了吧？我加上一次卷积再来一发倍增牛迭还是[tex]O(n \log n)[/tex]的。。
仿佛是在刻意逗我笑。
贴上自己丑的一比的代码：
void fft(int *x,int n,int fl=1){
for (int i=(n>>1),j=1;j<n;j++){
if (i<j) swap(x[i],x[j]); int k=(n>>1);
for (;k&i;i^=k,k>>=1); i^=k;
}
int now=0;
for (int m=2;m<=n;m<<=1){
int w; now++; if (fl==1) w=G[now]; else w=nG[now];
for (int i=0;i<n;i+=m){
int cur=1;
for (int j=i;j<i+(m>>1);j++){
int u=x[j],v=1ll*x[j+(m>>1)]*cur%mo;
x[j]=(u+v)%mo; x[j+(m>>1)]=(u-v+mo)%mo;
cur=1ll*cur*w%mo;
}
}
}
}
void getrev(int n){
if (n==1){
rev[0]=quick(fir[0],mo-2); return;
}
getrev((n+1)>>1); int len=1;
while ((1<<len)<=(n<<1)) len++;
for (int i=0;i<(1<<len);i++) A[i]=0;
for (int i=(n+1)>>1;i<(1<<len);i++) rev[i]=0;
for (int i=0;i<n;i++) A[i]=fir[i];
fft(A,1<<len); fft(rev,1<<len);
for (int i=0;i<(1<<len);i++) A[i]=(2-1ll*A[i]*rev[i]%mo+mo)%mo;
for (int i=0;i<(1<<len);i++) rev[i]=1ll*rev[i]*A[i]%mo; fft(rev,1<<len,-1);
for (int i=0;i<n;i++) rev[i]=1ll*rev[i]*two[len]%mo;
for (int i=n;i<(1<<len);i++) rev[i]=0;
}
void getln(int n){
getrev(n); for (int i=0;i<n-1;i++) rln[i]=1ll*fir[i+1]*(i+1)%mo;
int len=1; while ((1<<len)<=(n<<1)) len++;
for (int i=n-1;i<(1<<len);i++) rln[i]=0;
for (int i=n;i<(1<<len);i++) rev[i]=0;
fft(rln,1<<len); fft(rev,1<<len);
for (int i=0;i<(1<<len);i++) rln[i]=1ll*rln[i]*rev[i]%mo;
fft(rln,1<<len,-1);
for (int i=0;i<n;i++) rln[i]=1ll*two[len]*rln[i]%mo;
for (int i=n-1;i;i--) rln[i]=1ll*rln[i-1]*revnum[i]%mo;
rln[0]=0; for (int i=n;i<(1<<len);i++) rln[i]=0;
}
void getexp(int n){
if (n==1){
ans[0]=1; return;
}
getexp((n+1)>>1); int len=1;
while ((1<<len)<=(n<<1)) len++;
for (int i=0;i<n;i++) fir[i]=ans[i]; getln(n);
for (int i=n;i<(1<<len);i++) rln[i]=0;
for (int i=n;i<(1<<len);i++) fir[i]=0;
for (int i=0;i<(1<<len);i++) A[i]=0;
for (int i=0;i<n;i++) A[i]=C[i];
fft(A,1<<len); fft(fir,1<<len); fft(rln,1<<len);
for (int i=0;i<(1<<len);i++) ans[i]=(1ll*fir[i]*(1-rln[i]+A[i])%mo+mo)%mo;
fft(ans,1<<len,-1);
for (int i=n;i<(1<<len);i++) ans[i]=0;
for (int i=0;i<n;i++) ans[i]=1ll*ans[i]*two[len]%mo;
}
评论 (0)
