BZOJ 3527[Zjoi2014]力 FFT

題目連結:BZOJ3527

第一次學會如何寫數學公式，雖然隻是簡單的入門，但還是有點激動。。。

首先這個題很明顯是多項式乘法，但是強迫症的我過于糾結下标，以至于困惑了好久，簡直SB。

注:下标均從0開始。

現在進入正題。

Fj=∑i<jqiqj(i−j)2−∑i>jqiqj(i−j)2

現在我們轉化一下:記住下标從0開始

Fj=∑i=0j−1qiqj(i−j)2−∑i=j+1n−1qiqj(i−j)2

Ej=Fjqj=∑i=0j−1qi(i−j)2−∑i=j+1n−1qi(i−j)2

記 f(i)=qi ， g(i)=1i2 ，其中 g(0)=0 有

Ej=∑i=0j−1f(i)∗g(j−i)−∑i=j+1n−1f(i)∗g(j−i)

等号右邊第一個式子 ∑j−1i=0f(i)∗g(j−i)=∑ji=0f(i)∗g(j−i) ，就是卷積的形式，直接FFT算即可。

等号右邊第二個式子 ∑n−1i=j+1f(i)∗g(j−i)=∑n−1i=jf(i)∗g(j−i)=∑n−j−1i=0f(i+j)∗g(i)

倒序處理，令 h(n−1−i−j)=f(i+j) ，有 ∑n−j−1i=0f(i+j)∗g(i)=∑n−j−1i=0h(n−1−i−j)∗g(i) ，記為 Xi=∑n−j−1i=0h(n−1−i−j)∗g(i) ，則 Xn−j−1=∑ji=0h(j−i)∗g(i) ，也是卷積的形式，直接FFT即可。

最後 Ej=∑ji=0f(i)∗g(j−i)−Xn−j−1 ，其中 ∑ji=0f(i)∗g(j−i) 與 Xn−j−1 均可由FFT直接算的，題目即可解決。

代碼如下:

/**************************************************************
    Problem: 3527
    User: sfailsthy
    Language: C++
    Result: Accepted
    Time:3764 ms
    Memory:17680 kb
****************************************************************/

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
const int maxn;
const double PI =acos);

struct Complex{
    double real,image;
    Complex(){
        real=image;
    }

    Complex(double a,double b){
        real=a;
        image=b;
    }

    Complex operator + (const Complex &s) const{
        return Complex(real+s.real,image+s.image);
    }

    Complex operator - (const Complex &s) const{
        return Complex(real-s.real,image-s.image);
    }

    Complex operator * (const Complex &s) const{
        return Complex(real*s.real-image*s.image,real*s.image+image*s.real);
    }
}x1[maxn],x2[maxn],A[maxn];

int n,rev[maxn];
double a[maxn],b[maxn],c[maxn];

void init(int &len,int len1,int len2){
    int k,L,r,t;
    while(k*len1||k*len2){
        k<<;
        L++;
    }

    len=k;
    for(int i;i<len;i++){
        t=i,r,k=L;
        while(k--){
            r<<;
            r|=t;
            t>>;
        }
        rev[i]=r;
    }
}

void FFT(Complex x[],int len,int op){
    Complex u,t;
    for(int i;i<len;i++){
        A[rev[i]]=x[i];
    }
    for(int i;i<len;i++){
        x[i]=A[i];
    }

    for(int k;k<=len;k<<){
        Complex wn(cos*PI/k*op),sin*PI/k*op));
        for(int i;i<len;i+=k){
            Complex w);
            for(int j;j<k;j++){
                u=x[i+j];
                t=w*x[i+j+k];
                x[i+j]=u+t;
                x[i+j+k]=u-t;
                w=w*wn;
            }
        }
    }

    if(op=){
        for(int i;i<len;i++){
            x[i].real/=len;
        }
    }
}

int main(){
    scanf("%d",&n);
    for(int i;i<n;i++){
        scanf("%lf",&a[i]);
    }
    for(int i;i<n;i++){
        b[i]/i/i;
    }

    int len,len1=n,len2=n;
    init(len,len1,len2);
    for(int i;i<len;i++){
        x1[i]=Complex();
        x2[i]=Complex();
    }
    for(int i;i<len;i++){
       x1[i]=Complex(a[i]);
       x2[i]=Complex(b[i]);
    }

    FFT(x1,len);
    FFT(x2,len);
    for(int i;i<len;i++){
        x1[i]=x1[i]*x2[i];
    }
    FFT(x1,len);
    for(int i;i<n;i++){
        c[i]=x1[i].real;
    }

    for(int i;i<len;i++){
        x1[i]=Complex();
        x2[i]=Complex();
    }
    for(int i;i<n;i++){
        x1[i]=Complex(a[n-i]);
        if(i) x2[i]=Complex(b[i]);
    }
    FFT(x1,len);
    FFT(x2,len);
    for(int i;i<len;i++){
        x1[i]=x1[i]*x2[i];
    }
    FFT(x1,len);
    for(int i;i<n;i++){
        c[i]-=x1[n-i].real;
    }

    for(int i;i<n;i++){
        printf("%.3lf\n",c[i]);
    }
    return;
}