洛谷 P3338 [ZJOI2014]力（bzoj P3527 [ZJOI2014]力）

传送门

很难受，不会输公式，于是手写了一份，比较丑，希望不要介意。fft的套路果然深。。照着网上模板学习了一下。。

代码：

#include<iostream>
#include<stdio.h>
#include<string.h>
#include<math.h> 
using namespace std;
const int Maxn=1000005;
int n,N,L;
int rev[Maxn];
int dig[Maxn];
double p[Maxn];
struct complex
{
	double real,image;
	complex()
	{
		real=image=0.0;
	}
	complex(double a,double b)
	{
		real=a;
		image=b;
	}
}a[Maxn],b[Maxn],c[Maxn],A[Maxn],x,y,ans[Maxn];
complex operator +(complex a,complex b)
{
	return complex(a.real+b.real,a.image+b.image);
}
complex operator -(complex a,complex b)
{
	return complex(a.real-b.real,a.image-b.image);
}
complex operator *(complex a,complex b)
{
	return complex(a.real*b.real-a.image*b.image,a.real*b.image+a.image*b.real);
}
void FFT(complex a[],int flag)
{  
    for(int i=0;i<N;i++)
		A[i]=a[rev[i]];  
    for(int i=0;i<N;i++)
		a[i]=A[i];  
    for(int i=2;i<=N;i<<=1)
	{  
        complex wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));  
        for(int k=0;k<N;k+=i)
		{  
            complex w(1,0);  
            for(int j=0;j<i/2;j++)
			{
                x=a[k+j];  
                y=w*a[k+j+i/2];  
                a[k+j]=x+y;  
                a[k+j+i/2]=x-y;  
                w=w*wn;    
            }  
        }  
    }  
    if(flag==-1)for(int i=0;i<N;i++)a[i].real/=N;  
}  
int main()
{
	scanf("%d",&n);
	for(int i=0;i<n;i++)
		scanf("%lf",&p[i]);
	for(L=0,N=1;N<n;N<<=1,L++);
	N<<=1;
	L++;
	for(int i=0;i<N;i++)
	{
        int len=0;  
        for(int t=i;t;t>>=1)
			dig[len++]=t&1;  
        for(int j=0;j<L;j++)
			rev[i]=rev[i]*2+dig[j];  
    }
    for(int i=0;i<n;i++)
		a[i]=complex(p[i],0);
    for(int i=1;i<n;i++)
		b[i]=complex(1.0/i/i,0);
    FFT(a,1);
	FFT(b,1);
	for(int i=0;i<N;i++)
		c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)
		ans[i].real=c[i].real;
	memset(a,0,sizeof(a));  
    memset(b,0,sizeof(b));  
	for(int i=0;i<n;i++)
		a[i]=complex(p[n-i-1],0);
    for(int i=1;i<n;i++)
		b[i]=complex(1.0/i/i,0);
	FFT(a,1);
	FFT(b,1);
	for(int i=0;i<N;i++)
		c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<n;i++)
		ans[i].real-=c[n-i-1].real;
	for(int i=0;i<n;i++)
		printf("%f\n",ans[i].real);
	return 0;
}