• 【BZOJ4827】礼物(AHOI&HNOI2017)-FFT


    测试地址:礼物
    做法:本题需要用到FFT。
    为了方便讨论,我们把装饰物的编号设为0~n1,并把所有的调整操作都仅视为对第一串手环进行,显然这样并不失一般性。
    不难看出,第一串手环旋转k位,并增加c的亮度(如果原操作是对第二串手环亮度增加c,那么这里就是增加c)的差异值为:
    i=0n1(x(i+k)%n+cyi)2
    j=(i+k)%n,展开上述式子,最终得到差异值的式子为:
    (i=0n1xi2)+(i=0n1yi2)+2c(i=0n1xi)2c(i=0n1yi)+(nc2)2(i=0n1xjyi)
    我们发现差异值的最小值中c的取值一定在区间[m,m]内,而m只有100,所以考虑枚举c,那么上式中只有一个值取决于k
    i=0n1xjyi
    因为前面2的系数,所以我们要求这个东西的最大值,可是暴力求是O(n2)的,怎么办呢?
    实际上,答案就是两个向量A(Ai=xi)B(Bi=yn(i%n)1)(实际上就是yn1,yn2,...,y0的循环)的卷积C中的Cn1~C2n2这些位中的最大值。
    我们怎么能得到这个东西?因为上述从向量B到向量C中指定的那些位的变换,实际上就相当于右乘一个反循环矩阵X,即第一行为x0,x1,...,xn1,接下来每一行元素都较上一行左移了一位的矩阵,之所以称之为反循环矩阵是因为正的循环矩阵应该是每一行元素较上一行右移一位。因为我们不是很关心C中值的顺序,只需要以某种映射一一对应即可,所以我们把将矩阵第二行以下上下颠倒,变为正循环矩阵X,这样不会破坏一一对应的关系。通过某种规律,我们知道向量乘正循环矩阵可以用FFT优化,方法如下:
    将该正循环矩阵的第一行取出来,左右颠倒,然后向后循环,形成一个向量B(下标从0开始),然后将向量A(下标从0开始)和向量B做卷积,得到的向量C中的Cn1~C2n2就是所求答案。
    严谨的证明不作介绍,貌似也有前人研究过这种特殊的矩阵乘法,好像叫傅里叶对角化什么的,总之这样就把O(n2)的乘法优化成O(nlogn)了。而上面的答案很明显就可以看出是用这个方法推出的,所以是正确的。
    所以我们用O(nlogn)时间处理出i=0n1xjyi的最大值,用O(n)时间处理出xi,yi的和及平方和,然后O(m)枚举c来求最小差异值,总的复杂度和O(nlogn)同阶,完美解决了这一问题。
    以下是本人代码:

    #include <bits/stdc++.h>
    #define ll long long
    const double eps=0.5;
    const double pi=acos(-1.0);
    const ll inf=1000000000;
    using namespace std;
    ll n,m,suma=0,sumb=0,sumsq=0,r[500010];
    struct Complex
    {
        double x,y;
    }a[500010],b[500010];
    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 n,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("%lld%lld",&n,&m);
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        for(ll i=0;i<n;i++)
        {
            scanf("%lf",&a[i].x);
            suma+=(int)(a[i].x+eps);
            sumsq+=(ll)(a[i].x*a[i].x+eps);
        }
        for(ll i=0;i<n;i++)
        {
            scanf("%lf",&b[i].x);
            sumb+=(ll)(b[i].x+eps);
            sumsq+=(ll)(b[i].x*b[i].x+eps);
        }
        for(ll i=0;i<n-i-1;i++)
            swap(b[i].x,b[n-i-1].x);
        for(ll i=1;i<n;i++)
            b[n+i-1].x=b[i-1].x;
    
        ll bit=0,x=1;
        while (x<(n<<2)) bit++,x<<=1;
        r[0]=0;
        for(ll i=1;i<x;i++)
            r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
    
        FFT(a,(int)x,1);FFT(b,(int)x,1);
        for(ll i=0;i<x;i++) a[i]=a[i]*b[i];
        FFT(a,(int)x,-1);
    
        ll mx=0,ans=inf*inf;
        for(ll i=n-1;i<(n<<1)-1;i++)
            mx=max(mx,(ll)(a[i].x+eps));
        for(ll i=-m;i<=m;i++)
            ans=min(ans,sumsq+2*i*suma+n*i*i-2*mx-2*i*sumb);
        printf("%lld",ans);
    
        return 0;
    }
  • 相关阅读:
    xml 转义特殊字符
    MySql按周,按月,按日分组统计数据
    js jquery 动态添加表格
    详解Oracle partition分区表
    DBMS_NETWORK_ACL_ADMIN (OCP 053 第七题)
    RMAN 的优缺点及RMAN 备份及恢复步骤
    完美解决xhost +报错: unable to open display "" 装oracle的时候总是在弹出安装界面的时候出错
    Linux系统(X64)7 安装Oracle11g完整安装图文教程另附基本操作
    11g新特性-SQL Plan Management
    Oracle中事务处理控制用法
  • 原文地址:https://www.cnblogs.com/Maxwei-wzj/p/9793533.html
Copyright © 2020-2023  润新知