洛谷3338 [ZJOI2014]力

Luogu

BZOJ

分析

把条件展开,得到

$$ F_i=\sum\limits_{j=1}^{i-1}\frac{q_iq_j}{(i-j)^2}-\sum\limits_{j=i+1}^n\frac{q_iq_j}{(i-j)^2} $$

进一步得到

$$ E_i=\frac{F_i}{q_i}=\sum\limits_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum\limits_{j=i+1}^n\frac{q_j}{(i-j)^2} $$

设 $x_i=\frac{1}{i^2}$,得到

$$ E_i=\sum\limits_{j=1}^{i-1}q_jx_{i-j}-\sum\limits_{j=i+1}^nq_jx_{i-j} $$

设 $p_i=q_{n-i+1}$ ,得到

$$ E_i=\sum\limits_{j=1}^{i-1}q_jx_{i-j}-\sum\limits_{j=i+1}^np_{n-j+1}x_{i-j} $$

两个和式都是卷积的形式,FFT 求即可。

最后得到的结果再作差就是 $E_i$ 。

代码

代码中的f是上面的qg是上面的ph是上面的x

//It is made by M_sea
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;

const int MAXN=2000010;
const double PI=acos(-1.0);

struct Complex {
    double r,i;
    Complex(double x=0,double y=0) { this->r=x,this->i=y; }
    Complex operator +(const Complex x) { return Complex(r+x.r,i+x.i); }
    Complex operator -(const Complex x) { return Complex(r-x.r,i-x.i); }
    Complex operator *(const Complex x) { return Complex(r*x.r-i*x.i,r*x.i+i*x.r); }
} f[MAXN],g[MAXN],h[MAXN];

int n,limit;
int r[MAXN];

inline void FFT(Complex* a,int type) {
    for (re int i=0;i<limit;++i)
        if (i<r[i]) swap(a[i],a[r[i]]);
    for (re int i=1;i<limit;i<<=1) {
        Complex rot(cos(PI/i),type*sin(PI/i));
        for (re int j=0;j<limit;j+=i<<1) {
            Complex w(1,0);
            for (re int k=0;k<i;++k,w=w*rot) {
                Complex x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
}

int main() {
    scanf("%d",&n);
    for (re int i=1;i<=n;++i) scanf("%lf",&f[i].r);
    for (re int i=1;i<=n;++i) g[n-i+1].r=f[i].r;
    for (re int i=1;i<=n;++i) h[i].r=1.0/i/i;
    int l=0;
    for (limit=1;limit<=(n<<1);limit<<=1,++l);
    for (re int i=0;i<limit;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(f,1),FFT(g,1),FFT(h,1);
    for (re int i=0;i<=limit;++i) f[i]=f[i]*h[i],g[i]=g[i]*h[i];
    FFT(f,-1),FFT(g,-1);
    for (re int i=0;i<=limit;++i) f[i].r/=limit,g[i].r/=limit;
    for (re int i=1;i<=n;++i) printf("%.3lf\n",f[i].r-g[n-i+1].r);
    return 0;
}
最后修改:2019 年 09 月 24 日 08 : 19 PM

发表评论