ZJOI2014 力

一道简单的\(FFT\)题
题目链接

题意简述

给定一个公式\(E_i=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}\)
求\(E\).

---

解析

先把公式抄下来\(E_i=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}\)
我们令
\(A_i=q_i\)
\(i<0,B_i=-\frac{1}{i^2}\)
\(i=0,B_i=0\)
\(i>0,B_i=\frac{1}{i^2}\)
于是发现\(E_i=\sum_{j=0}^{n-1}A_j{B_{i-j}}\)
然后就是\(FFT\)板子题了.
注意读入的时候要平移一下,把\(B_i\)的\(i\)变成自然数.

代码如下

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#define N (600010)
#define inf (0x7f7f7f7f)
#define rg register int
#define Label puts("NAIVE")
#define spa print(' ')
#define ent print('\n')
#define rand() (((rand())<<(15))^(rand()))
typedef long double ld;
typedef long long LL;
typedef unsigned long long ull;
using namespace std;
const ld PI=3.14159265359;
struct com{
    ld a,b;
    com(){a=b=0;}
    com operator +(com x){
        com res;
        res.a=a+x.a,res.b=b+x.b;
        return res;
    }
    com operator -(com x){
        com res;
        res.a=a-x.a,res.b=b-x.b;
        return res;
    }
    com operator *(com x){
        com res;
        res.a=a*x.a-b*x.b,res.b=a*x.b+b*x.a;
        return res;
    }
}a[N],b[N];
int n,rev[N],Lim,len;
void FFT(com *a,int tp){
    for(int i=0;i<Lim;i++)
    if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=1;i<Lim;i<<=1){
        com w; w.a=cos(PI/i),w.b=tp*sin(PI/i);
        for(int R=i<<1,j=0;j<Lim;j+=R){
            com p; p.a=1;
            for(int k=j;k<j+i;k++,p=p*w){
                com x=a[k],y=p*a[k+i];
                a[k]=x+y,a[k+i]=x-y;
            }
        }
    }
    if(tp==-1)
    for(int i=0;i<Lim;i++)a[i].a=a[i].a/(ld)Lim;
}
int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    scanf("%Lf",&a[i].a);
    for(int i=-n+1;i<=n-1;i++){
        if(i==0)b[i+n-1].a=0;
        else b[i+n-1].a=1.0/(ld)i/(ld)i*((i<0)?(-1):1);
    }
    for(Lim=1;Lim<=(n*3);Lim<<=1)len++;
    for(int i=0;i<Lim;i++)
    rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    FFT(a,1),FFT(b,1);
    for(int i=0;i<Lim;i++)
    a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=n-1;i<=2*n-2;i++)
    printf("%.10Lf\n",a[i].a);
}

转载于:https://www.cnblogs.com/Romeolong/p/10075103.html