bzoj 3527 [快速傅里叶变换]

题目

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

$$F_j=\sum_{i< j}\frac{q_i\cdot q_j}{(i-j)^2}- \sum_{i> j}\frac{q_i\cdot q_j}{(i-j)^2}$$
令Ei=Fi/qi，求Ei.

输入

第一行一个整数n。
接下来n行每行输入一个数，第i行表示qi。
n≤100000，0 < qi < 1000000000

输出

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

样例输入

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

样例输出

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

分析

阅读题目,发现题目中给出的Ei=Fi/qi的条件可以化简,即:

$$E_j=\sum_{i< j}\frac{q_i}{(i-j)^2}- \sum_{i> j}\frac{q_i}{(i-j)^2}$$
这个式子很容易让我们将它与卷积联系起来,我们可以尝试构造两组多项式,使其乘积后系数中包含Ei,于是:
$$A(x)=q_0\cdot x_0+q_1\cdot x_1+q_2\cdot x_2+\dots +q_{n-1}\cdot x_{n-1} +0\cdot x_n + 0\cdot x_{n+1} + \dots 0\cdot x_{2*n-3}$$

$$B(x)=\frac{-1}{(n-1)^2}\cdot x_0 +\frac{-1}{(n-2)^2}\cdot x_1 + \dots+\frac{-1}{(-1)^2}\cdot x_{n-1}+0\cdot x_n+\frac{1}{1^2}\cdot x_{n+1}+\dots +\frac{1}{(n-2)^2}\cdot x_{2*n-3} +\frac{1}{(n-1)^2}\cdot x_{2*n-2}$$

这样令C(x)=A(x)$\times$B(x),C(x)的第n-1~2*n-2项的系数即为所求!
如果有读者不清楚FFT或者卷积的,可以参考我的博客:快速傅里叶变换(FFT)

完整代码

#include<bits/stdc++.h>
#define pi acos(-1.0)
#define maxn 500010
//#define DEBUG
using namespace std;
int n,m,l=0;
int rev[maxn];

struct Complex
{
    double real,imag;
    Complex() {real=0,imag=0;}
    Complex(double real,double imag):real(real),imag(imag) {}   
    Complex operator + (const Complex rhs) { return Complex(real+rhs.real,imag+rhs.imag); }
    Complex operator - (const Complex rhs) { return Complex(real-rhs.real,imag-rhs.imag); }
    Complex operator * (const Complex rhs) { return Complex(real*rhs.real-imag*rhs.imag,imag*rhs.real+real*rhs.imag); }
};
Complex a[maxn],b[maxn];

void pre()
{
    m=n*2;
    for(n=1;n<m;n<<=1) l++;
    for(int i=0;i<n;++i) rev[i]=(rev[i>>1]>>1|(1&i)<<(l-1));
#ifdef DEBUG
    cerr<<"n="<<n<<" m="<<m<<" l="<<l<<endl;
    for(int i=0;i<n;++i) cerr<<i<<"-->"<<rev[i]<<endl;
#endif
}

void FFT(Complex a[],int n,int flag)
{
    for(int i=0;i<n;++i)
        if(rev[i]<i) swap(a[rev[i]],a[i]);
    for(int i=2;i<=n;i<<=1)
    {
        Complex dw(cos(2*pi/i),sin(2*pi*flag/i));
        for(int j=0;j<n;j+=i)
        {
            Complex w(1,0);
            for(int k=0;k<(i>>1);k++,w=dw*w)
            {
                Complex u=a[j+k];
                Complex t=w*a[j+k+(i>>1)];
                a[j+k]=u+t;
                a[j+k+(i>>1)]=u-t; 
            }
        }
    }
    if(flag==-1) 
        for(int i=0;i<n;++i) a[i].real=a[i].real/n;
}

int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;++i)
    {
        double x;
        scanf("%lf",&x);
        a[i].real=x;
    }
    for(int i=0;i<n-1;++i) b[i].real=-1.0/(n-i-1)/(n-i-1);
    for(int i=n;i<2*n-1;++i) b[i].real=1.0/(i-n+1)/(i-n+1);
#ifdef DEBUG
    for(int i=0;i<2*n-1;++i) printf("%.4lf ",b[i].real);
#endif
    pre();
    FFT(a,n,1);
    FFT(b,n,1);
    for(int i=0;i<=n;++i) a[i]=a[i]*b[i];
    FFT(a,n,-1);
    for(int i=m/2-1;i<m-1;++i) printf("%lf\n",a[i].real);
    return 0;
}