bzoj3527: [Zjoi2014]力
题目链接
题解
大意：
\[F_j=\sum_{i<j}\frac{q_i q_j} {(i-j)^2}-\sum_{i>j}\frac{q_i q_j} {(i-j)^2}
\]
\[E_i = F_i / Q_i
\]
求E
化简F得
\[F_j=\sum_{i=0}^{j-1}\frac{q_i q_j}{(i-j)^2}-\sum_{i=j+1}^{n-1}\frac{q_i q_j}{(i-j)^2}
\]
带入\(F\),化简\(E\)式得到
\[E_j=\frac{F_j}{q_j}=\sum_{i=0}^{j-1}\frac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n-1}\frac{q_i}{(i-j)^2}
\]
令\(f(i)=q_i,g(i)=\frac{1}{i^2}\),得到
\[E_j=\sum_{i=0}^{j-1}f(i)*g(j-i)-\sum_{i=j+1}^{n-1}f(i)*g(j-i)
\]
对与\(\sum_{i=0}^{j-1}f(i)*g(j-i)\)卷积式，FFT即可
因为\(g(0) = 0\)
对于\(\sum_{i=j+1}^{n-1}f(i)*g(j-i)\)
\[\sum_{i=j+1}^{n-1}f(i)*g(j-i)=\sum_{i=j}^{n-1}f(i)*g(j-i)=\sum_{i=0}^{n-j-1}f(i+j)*g(i)
\]
令\(h(n-1-i-j)=f(i+j)\)
得到原式\(K_i = \sum_{i=0}^{n-j-1}h(n-1-i-j)*g(i)\)
那么\(K_{n-j-1}=\sum_{i=0}^{j}h(j-i)*g(i)\)
该式也化成了卷积式,FFT求解即可
那么
\[E_j=\sum_{i=0}^jf(i)*g(j-i)-K_{n-j-1}
\]
FFT求解
有点卡精度，处理初值g时要吧i*i转化为double，或者/i/i
PS,写FFT真的是BUG多多
代码
#include<cmath>
#include<cstdio>
#include<algorithm>
inline int read() {
int x = 0,f = 1;
char c = getchar();
while(c < '0' || c > '9') {if (c == '-') f = -1;c = getchar();}
while(c <= '9' && c >= '0') x = x * 10 + c - '0',c = getchar();
return x * f;
}
const int maxn = 300007;
const double pi = acos(-1.0);
double a[maxn];
struct Complex {
double x,y;
Complex (double x = 0,double y = 0) : x(x),y(y) {};
}f[maxn],g[maxn],c[maxn];
Complex operator + (Complex a,Complex b) {return Complex(a.x + b.x,a.y + b.y); };
Complex operator - (Complex a,Complex b) {return Complex(a.x - b.x,a.y - b.y); };
Complex operator * (Complex a,Complex b) {return Complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x); }
int n,m;
int l,r[maxn];
int limit = 1;
void FFT(Complex * A,int type) {
for(int i = 0;i < limit;++ i)
if(i < r[i]) std::swap(A[i],A[r[i]]);
for(int mid = 1;mid < limit;mid <<= 1) {
Complex wn (cos(pi / mid) , type * sin(pi / mid));
for(int R = mid << 1,j = 0;j < limit;j += R) {
Complex w(1,0);
for(int k = 0;k < mid;k ++,w = w * wn) {
Complex x = A[j + k],y = w * A[j + mid + k];
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
if(type == - 1) for(int i = 0;i < limit;++ i) A[i].x /= limit;
}
int main() {
n = read();double x;
//printf("%d\n",n);
for(int i = 0;i < n; ++ i)
scanf("%lf",&x),f[i].x = g[n - i - 1].x = x;
for(int i = 1; i < n;++ i)
c[i].x = 1.0 / (1.0* i * i);
while(limit <= n + n) limit <<= 1,l ++;
//printf("%d\n",limit);
for(int i = 0;i < limit;++ i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));//,printf("%d ",r[i]);;
//puts("");puts("");
FFT(f,1) ,FFT(g,1);FFT(c,1);
for(int i = 0;i < limit;++ i) f[i] = f[i] * c[i],g[i] = g[i] * c[i];
FFT(f,-1);FFT(g,-1);
for(int i = 0;i < n;++ i) printf("%.3lf\n",f[i].x - g[n - i -1].x);
return 0;
}
