BZOJ 3527: [Zjoi2014]力 [快速傅里叶变换]

n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

找到一个很详细的题解：http://blog.csdn.net/qq_33929112/article/details/54590319

花了两个小时来理解和写，1个小时查卡精度，又花了1个小时改模板卡洛谷的时限

以后都从0开始

观察这个式子

E=C-D

就是两个卷积啊....其中一个函数是1/i^2

规定1/0^2=0后，i=j也可以包括进来了然后前面那个是卷积标准形式，后面那个和上一题的技巧是一样的更改求和指标然后把q反转再反转D

注意：

傻逼题卡精度1/i/i才行

#include <iostream>
#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;
const int N=1e6+5;
const double PI=acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();}
    return x*f;
}
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
Vector conj(Vector a){return Vector(a.x,-a.y);}
 
struct FastFourierTransform{
    int n,rev[N];
    CD omega[N],omegaInv[N];
    void ini(int m){
        n=1;
        while(n<m) n<<=1;
        for(int k=0;k<n;k++) 
            omega[k]=CD(cos(2*PI/n*k),sin(2*PI/n*k)),
            omegaInv[k]=conj(omega[k]);
 
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
    }
    void transform(CD *a,CD *omega){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            for(CD *p=a;p!=a+n;p+=l)
                for(int k=0;k<m;k++){
                    CD t=omega[n/l*k]*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                }
        }
    }
    void DFT(CD *a,int flag){
        if(flag==1) transform(a,omega);
        else{
            transform(a,omegaInv);
            for(int i=0;i<n;i++) a[i].x/=(double)n;
        }
    }
}fft;
int n,m;
CD A1[N],A2[N],B[N];
double q;
int main(){
    //freopen("in","r",stdin);
    n=read();m=n+n-1;
    fft.ini(m);
    for(int i=0;i<n;i++) scanf("%lf",&q),A1[i].x=A2[n-i-1].x=q;
    fft.DFT(A1,1);fft.DFT(A2,1);
    for(int i=1;i<n;i++) B[i].x=1.0/i/i;
    fft.DFT(B,1);
     
    for(int i=0;i<fft.n;i++) A1[i]=A1[i]*B[i],A2[i]=A2[i]*B[i];
    fft.DFT(A1,-1);fft.DFT(A2,-1);
    for(int i=0;i<n;i++) printf("%.9lf\n",A1[i].x-A2[n-i-1].x);
}

#include <iostream>
#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;
const int N=1e6+5;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();}
    return x*f;
}
const double PI=acos(-1);
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

struct FastFourierTransform{
    int n,rev[N];
    void ini(int m){
        n=1;
        while(n<m) n<<=1;

        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
    }
    void DFT(CD *a,int flag){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            CD wn(cos(2*PI/l),flag*sin(2*PI/l));
            for(CD *p=a;p!=a+n;p+=l){
                CD w(1,0);
                for(int k=0;k<m;k++){
                    CD t=w*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                    w=w*wn;
                }
            }
        }
        if(flag==-1) for(int i=0;i<n;i++) a[i].x/=n;
    }
}fft;
int n,m;
CD A1[N],A2[N],B[N];
double q;
int main(){
    freopen("in","r",stdin);
    n=read();m=n+n-1;
    fft.ini(m);
    for(int i=0;i<n;i++) scanf("%lf",&q),A1[i].x=A2[n-i-1].x=q;
    fft.DFT(A1,1);fft.DFT(A2,1);
    for(int i=1;i<n;i++) B[i].x=1.0/i/i;
    fft.DFT(B,1);
    
    for(int i=0;i<fft.n;i++) A1[i]=A1[i]*B[i],A2[i]=A2[i]*B[i];
    fft.DFT(A1,-1);fft.DFT(A2,-1);
    for(int i=0;i<n;i++) printf("%lf\n",A1[i].x-A2[n-i-1].x);
}

BZOJ 3527: [Zjoi2014]力 [快速傅里叶变换]

3527: [Zjoi2014]力

Description

Input

Output