[洛谷P3338] [ZJOI2014]力

最新推荐文章于 2019-10-30 19:00:46 发布

转载最新推荐文章于 2019-10-30 19:00:46 发布 · 143 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/BCOI/p/10299474.html

本文解析洛谷P3338题目，通过数学变换将问题转化为卷积形式，利用FFT快速傅里叶变换求解电荷系统的电场强度。详细介绍了将原问题转换为多项式卷积的过程，并提供了C++代码实现。

洛谷题目链接:P3338 [ZJOI2014]力

题目描述

给出n个数qi，给出Fj的定义如下：

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

令Ei=Fi/qi，求Ei.

输入输出格式

输入格式：

第一行一个整数n。

接下来n行每行输入一个数，第i行表示qi。

输出格式：

n行，第i行输出Ei。

与标准答案误差不超过1e-2即可。

输入输出样例

输入样例#1：

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

输出样例#1：

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

说明

对于30%的数据，n≤1000。

对于50%的数据，n≤60000。

对于100%的数据，n≤100000，0<qi<1000000000。

[spj 0.01]

题解:

首先考虑化式子(话说这玩意求的好像是电场强度诶）.

\[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_j=\frac{F_j}{q_j}\]
\[E_j=\sum_{i=0}^{n-1}q_i*\frac{1}{(j-i)^2}-\sum_{i=j+1}^{n-1}q_i*\frac{1}{(i-j)^2}\]
不妨设\(b_i=\frac{1}{i^2}\),则有\[E_j=\sum_{i=0}^{n-1}q_i*b_{j-i}-\sum_{i=j+1}^{n-1}q_i*b_{i-j}\]
前面一半两个符号相乘的下标之和是一个常数,也就是如果将\(q,b\)看成多项式的话,那么这个乘积就可以做卷积.因为若\((a_0+a_1*x+a_2*x^2+...+a_{n-1}*x^{n-1})*(b_0+b_1*x+b_2*x^2+...+b_{n-1}*x^{n-1})=c_0+c_1*x+c_2*x^2+...+c^{2n-1}*x^{2n-1}\),则有\(c_k=\sum^k_{i=0}a_i*b_{k-i}\)也就是这样相乘可以使某一项的次数相同.

那么为了将后面的式子也转化成卷积的形式,我们可以将后面的\(q\)数组翻转一下,用\(q^{'}(q^{'}_i=q_{n-1-i})\)来代替,则有:\[E_j=\sum_{i=0}^{n-1}q_i*b_{j-i}-\sum^{n-1}_{i=j+1}q^{'}_{n-i-1}*b_{i-j}\]

那么后面那一半也变成了卷积的形式,就可以直接FFT求了.如果将\(b\)数组乘入\(q\)数组,那么最后的\(E_i\)对应着\(q_i-q^{'}_{n-i-1}\).

#include<bits/stdc++.h>
using namespace std;
const int N = 2e6+5;
const double eps = 1e-4;
const double pi = acos(-1.0);
typedef complex <double> comp;

int n, m, len = 0, r[N];
comp q1[N], q2[N], b1[N], b2[N];

void FFT(comp *A, int f){
    for(int i = 0; i < m; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(int i = 1; i < m; i <<= 1){
        for(int j = 0; j < m; j += (i<<1)){
            comp wi(cos(pi/i), sin(f*pi/i)), x, y, w(1, 0);
            for(int k = 0; k < i; k++, w *= wi){
                x = A[j+k], y = w*A[i+j+k];
                A[j+k] = x+y, A[i+j+k] = x-y;
            }
        }
    }
    if(f == -1) for(int i = 0; i < m; i++) A[i] /= m;
}

int main(){
    ios::sync_with_stdio(false);
    double x; cin >> n;
    for(int i = 0; i < n; i++) cin >> x, q2[n-i-1] = q1[i] = x;
    for(int i = 1; i <= n; i++) b1[i] = b2[i] = (double)(1.0/i/i);
    for(m = 1; m <= n*2; m <<= 1) len++;
    for(int i = 0; i < m; i++) r[i] = (r[i>>1]>>1) | ((i&1)<<(len-1));
    cerr << endl;
    FFT(q1, 1), FFT(b1, 1), FFT(q2, 1), FFT(b2, 1);
    for(int i = 0; i < m; i++) q1[i] *= b1[i], q2[i] *= b2[i];
    FFT(q1, -1), FFT(q2, -1);
    for(int i = 0; i < n; i++) printf("%.3lf\n", q1[i].real()-q2[n-i-1].real());
    return 0;
}

转载于:https://www.cnblogs.com/BCOI/p/10299474.html