• 【BZOJ3527】力(ZJOI2014)-FFT


    测试地址:
    做法:本题需要用到FFT。
    把题目所给式子中的qj除掉,我们发现题目要求:
    Ej=i<jqi(ij)2i>jqi(ij)2
    然而本人自己做到这里就并不会推了(即使在知道要使用FFT的情况下)。在本人看了题解之后,发现要根据一种新的思路来思考——向量卷积。
    根据向量卷积的定义,两个向量A,B的卷积C为:
    Ci=j+k=iAjbk
    为了方便讨论,下述的下标可能是负数。
    我们发现若令A,B为:
    Ai=qi,Bi=sign(i)×1i2
    其中sign(i)i的符号,即当i<0时,sign(i)=1,否则sign(i)=1
    那么不难看出Ei=Ci
    为了方便计算向量卷积,我们自然要把B的负数下标转化成非负下标,将下标均加上n1即可。于是用FFT计算向量卷积,然后取Cn~C2n1这些项即可,这些就是我们要求的E
    以下为本人代码(代码中A的下标从0开始,所以最后的向量中从第n-1项开始取):

    #include <bits/stdc++.h>
    using namespace std;
    int n,r[600010];
    struct Complex
    {
        double x,y;
    }a[600010],b[600010];
    const double pi=acos(-1.0);
    Complex operator + (Complex a,Complex b) {Complex s={a.x+b.x,a.y+b.y};return s;}
    Complex operator - (Complex a,Complex b) {Complex s={a.x-b.x,a.y-b.y};return s;}
    Complex operator * (Complex a,Complex b) {Complex s={a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};return s;}
    
    void FFT(Complex *a,int type)
    {
        for(int i=0;i<n;i++)
            if (i<r[i]) swap(a[i],a[r[i]]);
        for(int mid=1;mid<n;mid<<=1)
        {
            Complex W={cos(pi/mid),type*sin(pi/mid)};
            for(int l=0,r=mid<<1;l<n;l+=r)
            {
                Complex w={1.0,0.0};
                for(int k=0;k<mid;k++,w=w*W)
                {
                    Complex x=a[l+k],y=w*a[l+mid+k];
                    a[l+k]=x+y;
                    a[l+mid+k]=x-y;
                }
            }
        }
        if (type==-1)
        {
            for(int i=0;i<n;i++)
                a[i].x/=n;
        }
    }
    
    int main()
    {
        scanf("%d",&n);
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        for(int i=0;i<n;i++) scanf("%lf",&a[i].x);
        for(int i=0;i<(n<<1)-1;i++)
        {
            if (i<n-1) b[i].x=-1.0/((double)(i-n+1))/((double)(i-n+1));
            else if (i>n-1) b[i].x=1.0/((double)(i-n+1))/((double)(i-n+1));
                 else b[i].x=0.0;
        }
    
        int save,bit=0,x=1;
        while(x<=(n<<1)) x<<=1,bit++;
        r[0]=0;
        for(int i=1;i<x;i++)
            r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
        save=n;n=x;
    
        FFT(a,1),FFT(b,1);
        for(int i=0;i<n;i++) a[i]=a[i]*b[i];
        FFT(a,-1);
        for(int i=save-1;i<(save<<1)-1;i++)
            printf("%lf
    ",a[i].x);
    
        return 0;
    }
  • 相关阅读:
    Java EE
    C++基础学习(二)之判断
    C++基础学习(一)之循环
    AutoCAD完全卸载
    point-cloud-annotation-tool编译发烧记_UBUNTU
    Ubuntu18.04安装QT5
    Ubuntu安装VTK-8.1
    常用Latex公式代码表[持续更新]
    shp矢量文件批处理裁剪栅格影像_IDL/ENVI
    摄影测量基本原理
  • 原文地址:https://www.cnblogs.com/Maxwei-wzj/p/9793536.html
Copyright © 2020-2023  润新知