BZOJ3527 力【FFT】

最新推荐文章于 2019-06-01 15:28:03 发布

DenyTian

最新推荐文章于 2019-06-01 15:28:03 发布

阅读量200

点赞数 1

分类专栏： ====FFT

本文链接：https://blog.csdn.net/qq_34896694/article/details/85016170

版权

====FFT 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

题目链接：https://www.lydsy.com/JudgeOnline/problem.php?id=3527

题解：
易知 $E_j = \sum_{i<j}\frac{q_i}{(i-j)^2 }-\sum_{i>j} \frac{q_i}{(i-j)^2 }$
令 $p_i={1 \over i^2}$ 则原式变为 $\sum_{i=0}^{j-1}q_i*p_{i-j}+\sum_{i=j+1}^nq_j*p_{i-j}$ 注意到左边是一个卷积的形式（一开始对这里不太明白，应该是 $p_n*x^n+p_{n-1}*x^{n-1}+...+p_1*x+p_0)*(q_n*x^n+q_{n-1}*x^{n-1}+...+q_1*x+q_0)$ 乘起来就是）
对于右边我们可以进行一个变形，令 $i n v q [i] = q [n - i + 1]$ ，变成： $\sum_{i=0}^{j-1}invq_j*p_{i-j}$ （因为 $p$ 中有一个平方所以无所谓）
$F F T$ 两次即可

PS：注意到因为卷积的时候相当于最终结果是 $x$ 的 $2 n$ 次方级别的，所以预处理的时候要小心QAQ被坑了1h

// by Balloons
#include <cstdio>
#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
#define mpr make_pair
#define debug() puts("okkkkkkkk")
#define rep(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)

using namespace std;

typedef long long LL;

const int inf = 1 << 30;
const int maxn=4e5+5;
const double pi=acos(-1);

int n,m;
double p[maxn],pp[maxn],q[maxn];
int r[maxn];
double resa[maxn],resb[maxn];

struct complex{
    double x,y;
    complex(){x=0,y=0;}complex(double ax,double ay):x(ax),y(ay){}
}b[maxn<<1],c[maxn<<1];
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);}

void fft(complex *f,int op){
    for(int i=0;i<n;i++)
        if(i<r[i]){
            complex tt=f[i];
            f[i]=f[r[i]];
            f[r[i]]=tt;
        }
        
    for(int p=2;p<=n;p<<=1){
        int len=p/2;
        complex tmp(cos(pi/len),op*sin(pi/len));
        
        for(int k=0;k<n;k+=p){
            complex buf(1,0);
            for(int l=k;l<k+len;l++){
                complex tt=buf*f[l+len];
                f[l+len]=f[l]-tt;
                f[l]=f[l]+tt;
                buf=buf*tmp;
            }
        }
    }
}

complex tmpa[maxn],tmpb[maxn];
void work(double *a,double *b,double *res){
    for(int i=0;i<n;i++)tmpa[i]=complex(a[i],0),tmpb[i]=complex(b[i],0);
    fft(tmpa,1);fft(tmpb,1);
    for(int i=0;i<n;i++)tmpa[i]=tmpa[i]*tmpb[i];
    fft(tmpa,-1);
    for(int i=1;i<=m;i++)res[i]=(tmpa[i].x)/n;
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf",&p[i]);
    
    for(int i=1;i<=n;i++){
        q[i]=(double)(1.0/i/i);
        pp[i]=p[i];
    }
    reverse(pp+1,pp+n+1);
    
    for(m=n,n=1;n<=m;n<<=1);n<<=1;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
    
    work(p,q,resa);
    work(pp,q,resb);
    
    for(int i=1;i<=m;i++)printf("%.3f\n",resa[i]-resb[m-i+1]);

    return 0;
}