1. 程式人生 > >[bzoj2194]快速傅立葉之二_FFT

[bzoj2194]快速傅立葉之二_FFT

快速傅立葉之二 bzoj-2194

題目大意:給定兩個長度為$n$的序列$a$和$b$。求$c$序列,其中:$c_i=\sum\limits_{j=i}^{n-1} a_j\times b_{j-i}$。

註釋:$1\le n\le 10^5$,$0\le a_i,b_i\le 100$。


想法

顯然這是一道$FFT$裸題。

如圖:

上面的序列就是$a$序列,下面就是$b$序列。

左圖如題意,我們發現當把$b$序列翻轉之後就變成了右圖的樣子,我們設為$d$序列。

我們把$a$序列和$d$序列想象成兩個多項式,做多項式乘法之後,$c_i$就等於$\sum\limits_{j=i}^{n-1} a_j\times b_{n-j-1}$

這個卷積的形式我們可以通過$FFT$加速。

最後,附上醜陋的程式碼.... ...

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define N 100010 
using namespace std; typedef double db;
const db pi=acos(-1);
inline char nc() {static char *p1,*p2,buf[100000]; return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;}
int rd() {int x=0; char c=nc(); while(!isdigit(c)) c=nc(); while(isdigit(c)) x=(x<<3)+(x<<1)+(c^48),c=nc(); return x;}
struct cp
{
	db x,y;
	cp() {x=y=0;}
	cp(db x_,db y_) {x=x_,y=y_;}
	cp operator + (const cp &a) const {return cp(x+a.x,y+a.y);}
	cp operator - (const cp &a) const {return cp(x-a.x,y-a.y);}
	cp operator * (const cp &a) const {return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
}a[N<<2],b[N<<2];
void fft(cp *a,int len,int flg)
{
	int i,j,k,t;
	cp tmp,w,wn;
	for(i=k=0;i<len;i++)
	{
		if(i>k) swap(a[i],a[k]);
		for(j=len>>1;(k^=j)<j;j>>=1);
	}
	for(k=2;k<=len;k<<=1)
	{
		t=k>>1;
		wn=cp(cos(2*pi*flg/k),sin(2*pi*flg/k));
		for(i=0;i<len;i+=k)
		{
			w=cp(1,0);
			for(j=i;j<i+t;j++)
			{
				tmp=a[j+t]*w;
				a[j+t]=a[j]-tmp;
				a[j]=a[j]+tmp;
				w=w*wn;
			}
		}
	}
	if(flg==-1) for(i=0;i<len;i++) a[i].x/=len;
}
int main()
{
	int n=rd(); for(int i=0;i<n;i++) a[i].x=rd(),b[n-i-1].x=rd();
	int len=1; while(len<=(n<<1)) len<<=1;
	fft(a,len,1); fft(b,len,1);
	for(int i=0;i<len;i++) a[i]=a[i]*b[i];
	fft(a,len,-1);
	for(int i=0;i<n;i++) printf("%d\n",(int)(a[n+i-1].x+0.1));
	return 0;
}

小結:$FFT$真強...注意模板別背錯了。