bzoj 3527 [快速傅里叶变换]

题目

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

Fj=∑i<jqi⋅qj(i−j)2−∑i>jqi⋅qj(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的条件可以化简,即:

Ej=∑i<jqi(i−j)2−∑i>jqi(i−j)2 这个式子很容易让我们将它与卷积联系起来,我们可以尝试构造两组多项式,使其乘积后系数中包含Ei,于是: A(x)=q0⋅x0+q1⋅x1+q2⋅x2+⋯+qn−1⋅xn−1+0⋅xn+0⋅xn+1+…0⋅x2∗n−3

B(x)=−1(n−1)2⋅x0+−1(n−2)2⋅x1+⋯+−1(−1)2⋅xn−1+0⋅xn+112⋅xn+1+⋯+1(n−2)2⋅x2∗n−3+1(n−1)2⋅x2∗n−2

这样令C(x)=A(x) × 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; }