BZOJ3527: [Zjoi2014]力

给出n个数q_i，给出F_j的定义如下： 
 
令E_i=F_i/q_i，求E_i.

题解：

为什么要乘qi再除以qi。。。

我们发现i-j+j=i为定值，也就是这是个卷积。可以分开算，前半部分

令 a[i]=1/(i*i),b[i]=q[i],然后它们的卷积c[i]就是前半部分的值。

后半部分我们直接把b翻转一下再来一次就可以了。

注意不要给a添加无用的和不存在的项。

代码：比原来的pascal快了10s 233

#include<cstdio>

#include<cstdlib>

#include<cmath>

#include<cstring>

#include<algorithm>

#include<iostream>

#include<vector>

#include<map>

#include<set>

#include<queue>

#include<string>

#define inf 1000000000

#define maxn 500000+5

#define maxm 20000000+5

#define eps 1e-10

#define ll long long

#define pa pair<int,int>

#define for0(i,n) for(int i=0;i<=(n);i++)

#define for1(i,n) for(int i=1;i<=(n);i++)

#define for2(i,x,y) for(int i=(x);i<=(y);i++)

#define for3(i,x,y) for(int i=(x);i>=(y);i--)

#define mod 1000000007

using namespace std;

inline int read()

{

    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=10*x+ch-'0';ch=getchar();}

    return x*f;

}
struct cp
{
    double x,y;
    cp operator +(cp b){return (cp){x+b.x,y+b.y};}
    cp operator -(cp b){return (cp){x-b.x,y-b.y};}
    cp operator *(cp b){return (cp){x*b.x-y*b.y,x*b.y+y*b.x};}
};
double d[maxn],ans[maxn];
//const double PI=3.1415926535898;
const double PI=acos(-1.0);
cp a[maxn],b[maxn],c[maxn],y[maxn];
int n,nn,m,len,rev[maxn];
void fft(cp *x,int n,int flag)
{
    for0(i,n-1)y[rev[i]]=x[i];
    for0(i,n-1)x[i]=y[i];
    for(int m=2;m<=n;m<<=1)
    {
        cp wn=(cp){cos(2.0*PI/m*flag),sin(2.0*PI/m*flag)};
        for(int i=0;i<n;i+=m)
        {
            cp w=(cp){1,0};int mid=m>>1;
            for0(j,mid-1)
            {
                cp u=x[i+j],v=x[i+j+mid]*w;
                x[i+j]=u+v;x[i+j+mid]=u-v;
                w=w*wn;
            }
        }
    }
    if(flag==-1)for0(i,n-1)x[i].x/=n;
}

int main()

{

    freopen("input.txt","r",stdin);

    freopen("output.txt","w",stdout);

    nn=n=read();
    n=2*n-1;m=1;
    while(m<=n)m<<=1,len++;n=m;
    for0(i,nn-1)scanf("%lf",&d[i]);
    for0(i,n-1)a[i]=(cp){d[i],0};
    for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};
    for2(i,nn,n-1)b[i]=(cp){0,0};
    for0(i,n-1)
    {
        int x=i,y=0;
        for1(j,len)y<<=1,y|=x&1,x>>=1;
        rev[i]=y;
    }
    fft(a,n,1);fft(b,n,1);
    for0(i,n-1)c[i]=a[i]*b[i];
    fft(c,n,-1);
    for0(i,nn-1)ans[i]=c[i].x;
    for0(i,nn-1)a[i]=(cp){d[nn-1-i],0};
    for2(i,nn,n-1)a[i]=(cp){0,0};
    for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};
    for2(i,nn,n-1)b[i]=(cp){0,0};
    fft(a,n,1);fft(b,n,1);
    for0(i,n-1)c[i]=a[i]*b[i];
    fft(c,n,-1);
    for0(i,nn-1)ans[i]-=c[nn-1-i].x;
    for0(i,nn-1)printf("%.5f
",ans[i]);

    return 0;

}