【BZOJ3527】力（FFT）

【BZOJ3527】力（FFT）

题面

Description

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

\[Fj=\sum_{i<j}\frac{q_i q_j}{(i-j)^2 }-\sum_{i>j}\frac{q_i q_j}{(i-j)^2 } \]

令\(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

题解

首先就把\(qj\)直接除掉
于是式子变成了
\(\sum_{i<j} \frac{qi}{(i-j)^2}\)另一半同理
考虑这一半
qi分别要和\(\frac{1}{1^2},\frac{1}{2^2}\)等数相乘
于是做一遍FFT
系数分别为\(q_1,q_2....q_n\)
\(\frac{1}{1^2},\frac{1}{2^2},...\)
这要乘出来的系数和式子的前一半一一对应
然后把\(q\)反过来
变成\(q_n,q_{n-1}.....q_1\)
再做一遍FFT
就是式子后面的那一边
再一一对应减去就是答案

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
#include<complex>
using namespace std;
#define MAX 500000
const double Pi=acos(-1);
double Q[MAX],ans[MAX];
int n;
double sqr(double x){return x*x;}
int r[MAX];
complex<double> a[MAX],b[MAX];
int N,M,l;
void FFT(complex<double> *P,int opt)
{
    for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
    for(int i=1;i<N;i<<=1)
    {
        complex<double> W(cos(Pi/i),opt*sin(Pi/i));
        for(int p=i<<1,j=0;j<N;j+=p)
        {
            complex<double> w(1,0);
            for(int k=0;k<i;k++,w*=W)
            {
                complex<double> X=P[j+k],Y=w*P[j+k+i];
                P[j+k]=X+Y;P[j+k+i]=X-Y;
            }
        }
    }
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)scanf("%lf",&Q[i]);

	N=M=n-1;
	for(int i=0;i<=N;++i)a[i]=Q[i+1],b[i]=1.0/sqr(i+1);
	M+=N;
	for(N=1;N<=M;N<<=1)++l;
	for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	FFT(a,1);FFT(b,1);
	for(int i=0;i<=N;++i)a[i]=a[i]*b[i];
	FFT(a,-1);
	
	for(int i=2;i<=n;++i)ans[i]+=(double)(a[i-2].real()/N);

	for(int i=0;i<=N;++i)a[i].real()=b[i].real()=a[i].imag()=b[i].imag()=0;
	for(int i=0;i<n;++i)a[i]=Q[n-i],b[i]=1.0/sqr(i+1);

	FFT(a,1);FFT(b,1);
	for(int i=0;i<=N;++i)a[i]=a[i]*b[i];
	FFT(a,-1);

	for(int i=n-1;i;--i)ans[i]-=(double)(a[n-i-1].real()/N);

	for(int i=1;i<=n;++i)printf("%.3lf\n",ans[i]);
	return 0;
}