12
16
2014
0

论逗逼的自我修养之再探FFT

今天在北京某培训营中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;
}
Category: 论逗逼的自我修养系列 | Tags: | Read Count: 2349

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com