首页 > 代码库 > BZOJ 2194 快速傅立叶之二 ——FFT

BZOJ 2194 快速傅立叶之二 ——FFT

【题目分析】

    咦,这不是卷积裸题。

    敲敲敲,结果样例也没过。

    看看看,卧槽i和k怎么反了。

    艹艹艹,把B数组取个反。

    靠靠靠,怎么全是零。

    算算算,最终的取值范围算错了。

    交交交,总算是A掉了。

【代码】

#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>

#include <map>
#include <set>
#include <queue>
#include <string>
#include <iostream>
#include <algorithm>

using namespace std;

#define maxn 500005
#define inf 0x3f3f3f3f
#define ll long long 
#define cp Complex
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)

void Finout()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
//    freopen("out.txt","w",stdout);
    #endif
}

int Getint()
{
    int x=0,f=1; char ch=getchar();
    while (ch<‘0‘||ch>‘9‘) {if (ch==‘-‘) f=-1; ch=getchar();}
    while (ch>=‘0‘&&ch<=‘9‘) {x=x*10+ch-‘0‘; ch=getchar();}
    return x*f;
}

struct Complex{
	double x,y;
	Complex operator + (Complex a) { return (Complex) {x+a.x,y+a.y}; }
	Complex operator - (Complex a) { return (Complex) {x-a.x,y-a.y}; }
	Complex operator * (Complex a) { return (Complex) {x*a.x-y*a.y,x*a.y+y*a.x}; }
}a[maxn],b[maxn],c[maxn];

int n,m,len,rev[maxn],sum;
const double pi=acos(-1.0);

void FFT(Complex * x,int n,int f)
{
	F(i,0,n-1) if (rev[i]>i) swap(x[rev[i]],x[i]); //构造迭代的形式
	for (int m=2;m<=n;m<<=1)
	{
		Complex wn=(Complex){cos(2.0*pi/m*f),sin(2.0*pi/m*f)}; //当前的主单位根
		for (int i=0;i<n;i+=m)
		{
			Complex w=(Complex){1.0,0};//单位根 W0 
			for (int j=0;j<(m>>1);++j)
			{
				Complex u=x[i+j],v=x[i+j+(m>>1)]*w;
				x[i+j]=u+v; x[i+j+(m>>1)]=u-v;//蝴蝶变换,对称计算 
				w=w*wn;//下一个单位根 
			}
		}
	}
}

int main()
{
    Finout(); sum=n=Getint();
    F(i,0,n-1) scanf("%lf%lf",&a[i].x,&b[n-i-1].x);
	m=1; n=n*2-1;
	while (m<=n) m<<=1,len++; n=m;
	F(i,0,n-1)
	{
		int t=i,r=0;
		F(j,1,len) r<<=1,r|=t&1,t>>=1;
		rev[i]=r;
	}
	FFT(a,n,1); FFT(b,n,1);
	F(i,0,n-1) c[i]=a[i]*b[i];
	FFT(c,n,-1);
	F(i,0,n-1) (c[i].x/=n)+=0.4;
	F(i,sum-1,sum*2-2) printf("%lld\n",(ll)c[i].x);
}

  

BZOJ 2194 快速傅立叶之二 ——FFT