1. 程式人生 > >[BZOJ3527][Zjoi2014]力:FFT

[BZOJ3527][Zjoi2014]力:FFT

分析

整理得下式:

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

假設\(n=5\),考慮這兩個陣列:

\(a:q_1 \quad q_2 \quad q_3 \quad q_4 \quad q_5\)

\(b:-\frac{1}{16} \quad -\frac{1}{9} \quad -\frac{1}{4} \quad -\frac{1}{1} \quad 0 \quad \frac{1}{1} \quad \frac{1}{4} \quad \frac{1}{9} \quad \frac{1}{16}\)

容易發現\(E\)陣列是把\(a,b\)陣列看做多項式各項係數作卷積後一些項的係數。

FFT即可。

程式碼

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <complex>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,a) for(int i=head[(a)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
typedef std::complex<double> Complex;

const int MAXN=100005;
const double pi=3.14159265358979;
int n,m,len;
int rev[MAXN<<3];
Complex a[MAXN<<3],b[MAXN<<3];

inline void fft(Complex *c,int dft){
    rin(i,0,n-1) if(i<rev[i])
        std::swap(c[i],c[rev[i]]);
    for(int mid=1;mid<n;mid<<=1){
        int r=(mid<<1);
        Complex u=(Complex){cos(pi/mid),dft*sin(pi/mid)};
        for(int l=0;l<n;l+=r){
            Complex v=1;
            for(int i=0;i<mid;i++,v*=u){
                Complex x=c[l+i],y=v*c[l+mid+i];
                c[l+i]=x+y;
                c[l+mid+i]=x-y;
            }
        }
    }
    if(dft==-1) rin(i,0,n-1)
        c[i]/=n;
}

int main(){
    scanf("%d",&n);
    n--;
    rin(i,0,n){
        double x;
        scanf("%lf",&x);
        a[i]=x;
    }
    m=(n<<1);
    rin(i,0,m){
        if(i<n) b[i]=-1.0/(n-i)/(n-i);
        else if(i==n) b[i]=0;
        else b[i]=1.0/(i-n)/(i-n);
    }
    int nn=n;
    for(m+=n,n=1;n<=m;n<<=1) len++;
    rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|(i&1)<<(len-1));
    fft(a,1);
    fft(b,1);
    rin(i,0,n-1) a[i]*=b[i];
    fft(a,-1);
    n=nn;
    rin(i,n,n+n) printf("%.10lf\n",a[i].real());
    return 0;
}