BZOJ3527 力 【FFT】

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

题解:
易知 E j = ∑ i &lt; j q i ( i − j ) 2 − ∑ i &gt; j q i ( i − j ) 2 E_j = \sum_{i&lt;j}\frac{q_i}{(i-j)^2 }-\sum_{i&gt;j} \frac{q_i}{(i-j)^2 } Ej=i<j(ij)2qii>j(ij)2qi
p i = 1 i 2 p_i={1 \over i^2} pi=i21则原式变为 ∑ i = 0 j − 1 q i ∗ p i − j + ∑ i = j + 1 n q j ∗ p i − j \sum_{i=0}^{j-1}q_i*p_{i-j}+\sum_{i=j+1}^nq_j*p_{i-j} i=0j1qipij+i=j+1nqjpij注意到左边是一个卷积的形式(一开始对这里不太明白,应该是 ( 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 ) (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) (pnxn+pn1xn1+...+p1x+p0)(qnxn+qn1xn1+...+q1x+q0)乘起来就是)
对于右边我们可以进行一个变形,令 i n v q [ i ] = q [ n − i + 1 ] invq[i]=q[n-i+1] invq[i]=q[ni+1],变成: ∑ i = 0 j − 1 i n v q j ∗ p i − j \sum_{i=0}^{j-1}invq_j*p_{i-j} i=0j1invqjpij(因为 p p p中有一个平方所以无所谓)
F F T FFT FFT两次即可

PS:注意到因为卷积的时候相当于最终结果是 x x x 2 n 2n 2n次方级别的,所以预处理的时候要小心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;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值