Luogu P3338 [ZJOI2014]力

P3338 [ZJOI2014]力

题目大意

给出 \(n\) 个数 \(q_1,q_2, \dots q_n\)，定义

\[F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2} \]

\[E_i~=~\frac{F_i}{q_i} \]

对 \(1 \leq i \leq n\)，求 \(E_i\) 的值。

分析

我们来复现一下推导过程，熟悉一下。

\[E_j = \sum_{i<j}\frac{q_i}{(i-j)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2} \]

我们考虑转化，设\(f[i]=q_i,g[i]=\frac{1}{i^2}\)，所以原式变为

\[E_j=\sum_{i<j}f[i]*g[j-i]-\sum_{i>j}f[i]*g[i-j] \]

令f[0]=0,g[0]=0，考虑一下为什么这么定义？因为我们想让i=j的项无意义，这点对很多推导中的式子定义都有意义。

则我们可以改写式子

\[E_j=\sum_{i=0}^{j}f[i]*g[j-i]-\sum_{i=j}^{n}f[i]*g[i-j] \]

此时可以发现左边已经变为卷积形式了，考虑将右边转化为卷积形式，我们将右边展开

\[f[j]*g[0]+f[j+1]*g[1]+f[j+2]*g[2]+...+f[j+(n-j)]*g[n-j] \]

因此右边可以改写为

\[\sum_{i=0}^{n-j}f[j+i]*g[i] \]

接下来，我们用到一个技巧，翻转。我们引入f'[i]=f[n-i]，由此式子可以写为。

\[\sum_{i=0}^{n-j}f'[n-(j+i)]*g[i] \]

即

\[\sum_{i=0}^{n-j}f'[(n-j)-i)]*g[i] \]

令t=n-j，则式子的最终形式为

\[\sum_{i=0}^{t}f'[t-i]*g[i] \]

最后总结一下，式子即为

\[E_j = \sum_{i=0}^{j}f[i]*g[j-i] + \sum_{i=0}^{t}f[t-i]*g[i] \]

最后，我们分别求左边与右边的卷积，最后答案即为\(ans[j] = L_j-R_{n-j}\)

Ac_code

#include<bits/stdc++.h>
#define ios ios::sync_with_stdio(false); cin.tie(0), cout.tie(0)
using namespace std;
const int N = 3e5 + 10;
const double PI = acos(-1);
struct Complex
{
    double x,y;
    Complex(double x = 0,double y = 0) : x(x), y(y) {}
};

//复数乘法：模长相乘，幅度相加
Complex operator * (Complex J, Complex Q) {return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);}
Complex operator - (Complex J, Complex Q) {return Complex(J.x - Q.x, J.y - Q.y);}
Complex operator + (Complex J, Complex Q) {return Complex(J.x + Q.x, J.y + Q.y);}

namespace FFT
{
    typedef vector<Complex> poly;
    int R[N];//二进制翻转 二进制位数 补齐的2的整数幂N
    int FFT_init(int n)
    {
        int L = 0,limit = 1;
        while(limit<=n) limit <<= 1,L ++ ;
        // 补成2的整次幂，也就是N
        for(int i = 0; i < limit; ++ i)
            R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
        return limit;
    }
    void FFT(poly &A, int type,int limit)
    {
        A.resize(limit);
        for(int i = 0; i < limit; ++ i)
            if(i < R[i])
                swap(A[i], A[R[i]]);
            //i小于R[i]时才交换，防止同一个元素交换两次，回到它原来的位置。
        //从底层往上合并
        for(int mid = 1; mid < limit; mid <<= 1) {
            //待合并区间长度的一半，最开始是两个长度为1的序列合并,mid = 1;
            Complex wn(cos(PI / mid), type * sin(PI / mid));//单位根w_n^1;

            for(int len = mid << 1, pos = 0; pos < limit; pos += len) {
                //len是区间的长度，pos是当前的位置,也就是合并到了哪一位
                Complex w(1, 0);//幂,一直乘，得到平方，三次方...

                for(int k = 0; k < mid; ++ k, w = w * wn) {
                    //只扫左半部分，蝴蝶变换得到右半部分的答案,w 为 w_n^k
                    Complex x = A[pos + k];//左半部分
                    Complex y = w * A[pos + mid + k];//右半部分
                    A[pos + k] = x + y;//左边加
                    A[pos + mid + k] = x - y;//右边减
                }
            }
        }
        if(type == 1) return ;
        for(int i = 0; i < limit; ++ i)
            A[i].x /= limit;
            //最后要除以limit也就是补成了2的整数幂的那个N，将点值转换为系数
            //（前面推过了点值与系数之间相除是N）
    }
    poly poly_mul(poly A,poly B)
    {
        int deg = A.size() + B.size() - 1;
        int limit = FFT_init(deg);
        poly C(limit);
        FFT(A,1,limit),FFT(B,1,limit);
        for(int i=0;i<limit;++i) C[i] = A[i]*B[i];
        // cout<<cnt<<endl;
        FFT(C,-1,limit);
        C.resize(deg);
        return C;
    }
};

using FFT::poly;
using FFT::poly_mul;

int main()
{
    int n;scanf("%d",&n);
    poly a(n+1),b(n+1),c(n+1);
    for(int i=1;i<=n;i++) 
    {
        scanf("%lf",&a[i].x);
        c[n-i].x = a[i].x;
        b[i].x = 1.0/i/i;
    }
    a = poly_mul(a,b,1);
    c = poly_mul(b,c,2);
    for(int i=1;i<=n;i++) printf("%.3lf\n",a[i].x-c[n-i].x);
    return 0;
}