【BZOJ-3527】力 FFT

3527: [Zjoi2014]力

Time Limit: 30 Sec Memory Limit: 256 MBSec Special Judge
Submit: 1544 Solved: 899
[Submit][Status][Discuss]

Description

给出n个数qi，给出Fj的定义如下：

令Ei=Fi/qi，求Ei.

Input

第一行一个整数n。

接下来n行每行输入一个数，第i行表示qi。

n≤100000，0<qi<1000000000

Output

n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

HINT

Source

Solution

一道裸的FFT，我瞪了快一节课...

先两边同除$p_{j}$就可以直接得到$E_{j}$的关于$q$的关系..然后就可以看出是一个卷积的形式了..

主要是光想直接$A\bigotimes B=E$，实际上正反求一下相减就好..

然后这里discuss里提醒了我一件事..爆int

Code

#include<iostream>

#include<cstdio>

#include<algorithm>

#include<cmath>

#include<cstring>

using namespace std;

struct Complex{

	double r,i;

	Complex(double R=0.0,double I=0.0) {r=R; i=I;}

	Complex operator + (const Complex & A) const {return Complex(r+A.r,i+A.i);}

	Complex operator - (const Complex & A) const {return Complex(r-A.r,i-A.i);}

	Complex operator * (const Complex & A) const {return Complex(r*A.r-i*A.i,r*A.i+i*A.r);}

};

#define MAXN 600010

#define Pai acos(-1.0)

Complex A[MAXN],B[MAXN],C[MAXN],D[MAXN];

int N,len;

double a[MAXN];

inline void Prework()

{

	len=1;

	while (len<((N-1)<<1)) len<<=1;

	for (int i=0; i<=N-1; i++) A[i]=Complex(a[i],0);

	for (int i=N; i<len; i++) A[i]=Complex(0,0);

	for (int i=0; i<=N-1; i++) B[i]=Complex(a[N-1-i],0);

	for (int i=N; i<len; i++) B[i]=Complex(0,0);

	for (int i=0; i<=N-1; i++) if (i) D[i]=C[i]=Complex(1.0/i/i,0);

	for (int i=N; i<len; i++) D[i]=C[i]=Complex(0,0);

}

inline void Rader(Complex *x)

{

    for (int i=1,j=len>>1,k; i<len-1; i++)

        {

            if (i<j) swap(x[i],x[j]);

            k=len>>1;

            while (j>=k) j-=k,k>>=1;

            if (j<k) j+=k;

        }

}

inline void DFT(Complex *x,int opt)

{

    Rader(x);

    for (int h=2; h<=len; h<<=1)

        {

            Complex Wn( cos(opt*2*Pai/h),sin(opt*2*Pai/h) );

            for (int i=0; i<len; i+=h)

                {

                    Complex W(1,0);

                    for (int j=i; j<i+h/2; j++)

                        {

                            Complex u=x[j],t=W*x[j+h/2];

                            x[j]=u+t; x[j+h/2]=u-t;

                            W=W*Wn;

                        }

                }

        }

    if (opt==-1)

        for (int i=0; i<len; i++) x[i].r/=len;

}

inline void FFT(Complex *x,Complex *y)

{

	DFT(x,1); DFT(y,1);

	for (int i=0; i<len; i++) x[i]=x[i]*y[i];

	DFT(x,-1);

}

int main()

{

	scanf("%d",&N);

	for (int i=0; i<N; i++) scanf("%lf",&a[i]);

	Prework();

//	for (int i=0; i<len; i++) printf("%.6lf\n",A[i].r); puts("=================");

//	for (int i=0; i<len; i++) printf("%.6lf\n",B[i].r); puts("=================");

//	for (int i=0; i<len; i++) printf("%.6lf\n",C[i].r); puts("=================");

	FFT(A,C); FFT(B,D);

	for (int i=0; i<N; i++) printf("%.3lf\n",A[i].r-B[N-1-i].r);

	return 0;

}

秒客网