博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
「洛谷3338」「ZJOI2014」力【FFT】
阅读量:4581 次
发布时间:2019-06-09

本文共 2830 字,大约阅读时间需要 9 分钟。

题目链接

题解

首先我们需要对这个式子进行化简,否则对着这么大一坨东西只能暴力。。。


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

根据题目给出的定义,带入\(E\)
\[E_i=\sum_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^{n}\frac{q_j}{(j-i)^2}\]
形式稍微改变了一下,本质一样
需要处理这个分母上的两个不好看的小东西
\(f_i=\frac{1}{i^2}\)
那么可以得到原式变成
\[E_i=\sum^{i-1}_{j=1}q_jf_{i-j}-\sum^n_{j=i+1}q_jf_{j-i}\]
将这个式子拆成两部分
先讨论\(-\)号前面一部分
可以发现,减号前一部分就是一个卷积的形式,那么直接FFT套上去就可以了。
对于后一部分,认真观察可以发现
形式和前一部分非常相似,其实就是相反了一下。
我们将原来的数组倒序存储一下,和前一部分做一样的FFT就可以了。

代码

#include 
#define ms(a, b) memset(a, b, sizeof(a))#define ll long long#define ull unsigned long long#define ms(a, b) memset(a, b, sizeof(a))#define inf 0x3f3f3f3f#define db double#define Pi acos(-1)#define eps 1e-8#define N 400005using namespace std;template
T sqr(T x) { return x * x; }template
void read(T &x) { x = 0; T fl = 1; char ch = 0; for (; ch < '0' || ch > '9'; ch = getchar()) if (ch == '-') fl = -1; for (; ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48); x *= fl;}template
T power(T x, T y, T mod) { T res = 1; for (; y; y >>= 1) { if (y & 1) res = (res * x) % mod; x = (x * x) % mod; } return res; }template
void write(T x) { if (x < 0) x = -x, putchar('-'); if (x > 9) write(x / 10); putchar(x % 10 + '0'); }template
void writeln(T x) { write(x); puts(""); }struct Complex { db x, y; Complex (db xx = 0, db yy = 0) { x = xx; y = yy; } Complex operator + (Complex B) { return Complex(x + B.x, y + B.y); } Complex operator - (Complex B) { return Complex(x - B.x, y - B.y); } Complex operator * (Complex B) { return Complex(x * B.x - y * B.y, x * B.y + y * B.x); }}a[N], b[N], c[N], d[N];int r[N];int n, m, limit, l;void FFT(Complex *a, int type) { for (int i = 0; i < limit; i ++) if (i < r[i]) 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 /= 1.0 * limit;}int main() { read(n); for (int i = 1; i <= n; i ++) { db x; scanf("%lf", &x); a[i].x = c[n - i + 1].x = x; b[i].x = d[i].x = 1.0 / sqr(i * 1.0); } limit = 1; while (limit <= (n << 1)) limit <<= 1, l ++; for (int i = 0; i < limit; i ++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); FFT(a, 1); FFT(b, 1); for (int i = 0; i < limit; i ++) a[i] = a[i] * b[i]; FFT(a, -1); FFT(c, 1); FFT(d, 1); for (int i = 0; i < limit; i ++) c[i] = c[i] * d[i]; FFT(c, -1); for (int i = 1; i <= n; i ++) printf("%.3lf\n", a[i].x - c[n - i + 1].x); return 0;}

转载于:https://www.cnblogs.com/chhokmah/p/10785747.html

你可能感兴趣的文章
jmeter- 性能测试3:聚合报告(Aggregate Report )
查看>>
JavaScript高级程序设计---学习笔记(二)
查看>>
vim 插件的学习
查看>>
Uncaught SyntaxError: Unexpected token ILLEGAL
查看>>
一个预处理定义的问题
查看>>
ANDROID L——Material Design综合应用(Demo)
查看>>
自我介绍以及关于软件工程的问题
查看>>
struts (一)
查看>>
【新番推荐】工作细胞
查看>>
NYOJ 16 矩形嵌套
查看>>
Leetcode中的SQL题目练习(二)
查看>>
dubbo 集群容错源码
查看>>
Collection接口的子接口——Queue接口
查看>>
LINUX安装NGINX
查看>>
服务器启动项目抛错 没有到主机的路由
查看>>
python_85_sys模块
查看>>
第九周动手动脑
查看>>
HDU 1811 Rank of Tetris
查看>>
winform 获取当前名称
查看>>
opengl glut vs2013配置
查看>>