[ZJOI2014]力
题目大意:
给定(q_{1sim n}(nle10^5)),求(g_i=sum_{j<i}frac{q_j}{(i-j)^2}+sum_{i>j}frac{q_j}{(i-j)^2})。
思路:
令(f_i=frac1{i^2}),则
[g_i=sum_{j=0}^{i-1}q_if_{i-j}-sum_{j=i+1}^nq_jf_{j-i}
]
令(p)为(q)的翻转,(a_i)表示前半个和式,(b_i)表示后半个和式翻转,则
[a_i=sum_{j=0}^{i-1}q_if_{i-j}
]
[b_i=sum_{j=0}^{i-1}p_if_{i-j}
]
用FFT即可求出。
时间复杂度(mathcal O(nlog n))。
源代码:
#include<cmath>
#include<cstdio>
#include<cctype>
#include<complex>
#include<algorithm>
inline int getint() {
register char ch;
while(!isdigit(ch=getchar()));
register int x=ch^'0';
while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
return x;
}
inline double getdouble() {
double x;
scanf("%lf",&x);
return x;
}
const int N=262144;
const double pi=M_PI;
typedef std::complex<double> comp;
int lim;
double ans[N];
comp q[N],ru[N],iru[N],p[N],f[N],a[N],b[N];
inline void init_ru(const int &n) {
for(register int i=0;i<n;i++) {
ru[i]=(comp){cos(2*pi*i/n),sin(2*pi*i/n)};
iru[i]=conj(ru[i]);
}
}
inline void dft(comp f[],comp w[],const int &n) {
for(register int i=0,j=0;i<n;i++) {
if(i>j) std::swap(f[i],f[j]);
for(register int l=n>>1;(j^=l)<l;l>>=1);
}
for(register int i=2;i<=n;i<<=1) {
const int m=i>>1;
for(register int j=0;j<n;j+=i) {
for(register int k=0;k<m;k++) {
const comp z=f[j+m+k]*w[n/i*k];
f[j+m+k]=f[j+k]-z;
f[j+k]+=z;
}
}
}
}
int main() {
const int n=getint();
for(register int i=0;i<n;i++) {
q[i]=getdouble();
}
for(lim=1;lim<n;lim<<=1);lim<<=1;
init_ru(lim);
std::copy(&q[0],&q[lim],&p[0]);
std::reverse(&p[0],&p[lim]);
for(register int i=1;i<n;i++) {
f[i]=1./i/i;
}
dft(q,ru,lim);
dft(p,ru,lim);
dft(f,ru,lim);
for(register int i=0;i<lim;i++) {
a[i]=q[i]*f[i];
b[i]=p[i]*f[i];
}
dft(a,iru,lim);
dft(b,iru,lim);
for(register int i=0;i<lim;i++) {
a[i]/=lim;
b[i]/=lim;
}
std::reverse(&b[0],&b[lim]);
for(register int i=0;i<n;i++) {
printf("%f
",a[i].real()-b[i].real());
}
return 0;
}