bzoj3527: [Zjoi2014]力

Posted 2020-11-01 Nico&11101001

题目链接

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;
}