[ZJOI2014]力【FFT】

[ZJOI2014]力

已知\[E(j)=\sum_{i < j} \frac{q_i}{(i-j)^2}-\sum_{i > j} \frac{q_i}{(i-j)^2}\]， 求\(E\)

令\(a(i)=q_i ,b(i)=i^{-2}\)

对于减号左边，\[L(j)=\sum_{i < j}a(i)b(j-i)\]

对于减号右边，\[R(j)=-\sum_{i > j}a(i)b(i-j)\]

令\[c(i)=b(i-n)=(i-n)^{-2}, i > n\]， 则\[L(j)=\sum_{i < j}a(i)b(j-i+n-n)=\sum_{i < j} a(i)c(j-i+n)=(a*c)(n+j)\]

令\[c(i)=-b(n-i)=-(n-i)^{-2}, i < n\]，即\[c(n-i)=-b(i)=-i^{-2}\],则\[R(j)=\sum_{i > j}a(i)b(i-j)=\sum_{i > j}a(i)c(n-i+j)=(a*c)(n+j)\]

特殊的，令\(c(n)=0\)

则\[E(j)=(a*c)(n+j)\]

求出\(a\)和\(c\)的卷积就好了

const double pi=acos(-1.0);
const int Maxn=2e6+10;
int n, l, r[Maxn], lim=1;

struct complex{double a, b; complex(double a=0, double b=0):a(a), b(b){}} f[Maxn], g[Maxn];
complex operator +(complex &A, complex &B){return complex(A.a+B.a, A.b+B.b);}
complex operator -(complex &A, complex &B){return complex(A.a-B.a, A.b-B.b);}
complex operator *(complex &A, complex &B){return complex(A.a*B.a-A.b*B.b, A.b*B.a+A.a*B.b);}

void FFT(complex *A, int tp){
    for(int i=0; i < lim; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    complex wn, w, x, y;
    for(int mid=1; mid < lim; mid<<=1){ wn=complex(cos(pi/mid), tp*sin(pi/mid));
        for(int j=0, R=mid<<1; j < lim; j+=R){ w=complex(1, 0);
            for(int k=0; k < mid; k++, w=w*wn) x=A[j+k], y=w*A[j+k+mid], A[j+k]=x+y, A[j+k+mid]=x-y;
        }
    }
}
void solve(){
    n=read();
    for(int i=1; i <= n; i++) scanf("%lf", &f[i].a); g[n].a=0;
    for(int i=0; i < n; i++) g[i].a=-1.0/(n-i)/(n-i);
    for(int i=n+1; i <= n*2; i++) g[i].a=1.0/(i-n)/(i-n);
    while(lim <= n*3) lim <<= 1, l++;
    for(int i=0; i < lim; i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
    FFT(f, 1); FFT(g, 1); for(int i=0; i < lim; i++) f[i]=f[i]*g[i];
    FFT(f, -1); for(int i=n+1; i <= n+n; i++) printf("%lf\n", f[i].a/lim);
}

原文地址：https://www.cnblogs.com/zerolt/p/9297922.html