BZOJ3527 推出卷积公式FFT求值

传送门:https://www.lydsy.com/JudgeOnline/problem.php?id=3527

题意:

\(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_i=F_i/q_i\)

题解:

推公式:

\[E_i=F_i/q_i\\
E_i=\sum_{j=i}^{n}\frac{q_j}{(i-j)^2}-\sum_{j=0}^{i}\frac{q_j}{(i-j)^2}\\
设函数f(i)为q_i,g(i)为(i)^2\\
\sum_{i=0}^{j} f_{i} * g_{j-i}\\
\begin{array}{l}{\sum_{i>j} \frac{q_{i}}{(i-j)^{2}}=\sum_{i=j}^{n} \frac{q_{i}}{(i-j)^{2}}} {=\sum_{i=0}^{n-j} \frac{q_{n-i}}{(j-i)^{2}}=\sum_{i=0}^{n-j} f_{n-i} * g_{i-j}}\end{array}
\]

于是我们求出f和g 后fft,然后求值即可

代码:

#include <set>
#include <map>
#include <cmath>
#include <cstdio>
#include <string>
#include <vector>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long LL;
typedef pair<int, int> pii;
typedef unsigned long long uLL;
#define ls rt<<1
#define rs rt<<1|1
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
#define bug printf("*********\n")
#define FIN freopen("input.txt","r",stdin);
#define FON freopen("output.txt","w+",stdout);
#define IO ios::sync_with_stdio(false),cin.tie(0)
#define debug1(x) cout<<"["<<#x<<" "<<(x)<<"]\n"
#define debug2(x,y) cout<<"["<<#x<<" "<<(x)<<" "<<#y<<" "<<(y)<<"]\n"
#define debug3(x,y,z) cout<<"["<<#x<<" "<<(x)<<" "<<#y<<" "<<(y)<<" "<<#z<<" "<<z<<"]\n"
const int maxn = 1e6 + 5;
const int INF = 0x3f3f3f3f;
const int mod = 1e9 + 7;
const double Pi = acos(-1.0);
LL quick_pow(LL x, LL y) {
LL ans = 1;
while(y) {
if(y & 1) {
ans = ans * x % mod;
} x = x * x % mod;
y >>= 1;
} return ans;
}
struct complex {
double x, y;
complex(double xx = 0, double yy = 0) {
x = xx, y = yy;
}
} f[maxn], f1[maxn], g[maxn];
complex operator + (complex a, complex b) {
return complex(a.x + b.x, a.y + b.y);
}
complex operator - (complex a, complex b) {
return complex(a.x - b.x, a.y - b.y);
}
complex operator * (complex a, complex b) {
return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
} int n, m;
int l, r[maxn];
int limit = 1;
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 + k + mid] = x - y;
}
}
}
}
int ans[maxn];
char numA[maxn], numB[maxn];
int main() {
#ifndef ONLINE_JUDGE
FIN
#endif
int n;
while(scanf("%d", &n) != EOF) {
n--;
m = n * 2;
for(int i = 0; i <= n; i++) {
scanf("%lf", &f[i].x);
// debug1(f[i].x);
f1[n - i].x = f[i].x;
}
// bug;
for(int i = 1; i <= n; i++) {
g[i].x = (double)(1.0 / i / i);
}
while(limit <= m) limit <<= 1, l++;
for(int i = 0; i <= limit; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
fft(f, 1);
fft(f1, 1);
fft(g, 1);
for(int i = 0; i <= limit; i++) {
f[i] = f[i] * g[i];
}
for(int i = 0; i <= limit; i++) {
g[i] = f1[i] * g[i];
}
fft(f, -1);
fft(g, -1);
for(int i = 0; i <= limit; i++) {
f[i].x = f[i].x / limit;
}
for(int i = 0; i <= limit; i++) {
g[i].x = g[i].x / limit;
}
int t = m / 2;
for(int i = 0; i <= t; i++) {
printf("%.3f\n", f[i].x - g[n - i].x);
} }
return 0;
}

最新文章

  1. C++ TR1 Function Bind
  2. RobotFramework自动化测试之环境搭建安装教程(一)
  3. Linux Kernel CMPXCHG函数分析
  4. java:IO-读写大文件
  5. Android 常用UI控件之TabHost(2)简单示例
  6. OC加强-day06
  7. 在Ubuntu 12.10 上安装部署Openstack
  8. AppDomain.CurrentDomain.GetAssemblies()
  9. Scala学习笔记--Actor和并发
  10. Immutable 详解及 React 中实践
  11. 蓝桥杯-有理数类-java
  12. 关于对JavaScript待于完善的一些知识点
  13. ArcSDE 10.1安装、配置、连接
  14. ynoi2018
  15. groovy.lang.GroovyRuntimeException: Conflicting module versions
  16. maven 相关问题
  17. Final互评------《I do》---- 二次元梦之队
  18. 大数据 -- Spark
  19. R语言学习——欧拉计划(3)Largest prime factor 求最大质因数
  20. 详解MathType快捷键使用技巧

热门文章

  1. Kubernetes1.3新特性:rktnetes
  2. Android开发-API指南-&lt;activity-alias&gt;[原创译文]
  3. BZOJ 3884 上帝与集合的正确用法题解
  4. 18-2 djanjo中间件和orm多对多操作,以及ajax
  5. jmter对于函数的处理
  6. HZOJ 走格子
  7. 如何使用jmeter调用soap协议
  8. hdu 1384 Intervals (差分约束)
  9. Java 内存模型及GC原理 (转)
  10. C#中的?操作符