项链 [FFT]

题面

1226867-20180929110816435-619625003.png

思路

这题很像bzoj4827礼物

还是一样的思路,我们把$y$倍长,$y[i+k]=y[i]+n$

然后令$f(s,c)$表示从$y$的第$s$个开始匹配,位置偏移量为$c$的答案

可以得到$f(s,c)=\sum_{i=0}^{n-1}(x_i-y_{i+s}+c)^2=\sum_{i=0}^{n-1}(x_i^2+y_{i+s}+c^2+2x_ic-2y_{i+s}x-2x_iy_{i+s})$

我们可以把右边这个式子视为关于$c$的二次函数

用$FFT$可以快速得到不同的$s$下,$2x_iy_{i+s}$这一项的值,又因为其他系数都是确定的,所以我们可以用$O(k\log k)$确定不同的s对应的c的函数

然后可以算出来这个函数的最佳取值,再从所有最佳取值中选出最优解即可

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
inline int read(){
    int re=0,flag=1;char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-') flag=-1;
        ch=getchar();
    }
    while(isdigit(ch)) re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
    return re*flag;
}
namespace FFT{
    struct complex{
        double x,y;
        complex(double xx=0,double yy=0){x=xx;y=yy;}
        inline complex operator +(complex &b){return complex(x+b.x,y+b.y);}
        inline complex operator -(complex &b){return complex(x-b.x,y-b.y);}
        inline complex operator *(complex &b){return complex(x*b.x-y*b.y,x*b.y+y*b.x);}
    }A[100010],B[100010];
    int cnt,lim,r[100010];
    const double pi=acos(-1.0);
    void fft(complex *a,double type){
        int i,j,k,mid;complex x,y,w,wn;
        for(i=0;i<lim;i++) if(i<r[i]) swap(a[i],a[r[i]]);
        for(mid=1;mid<lim;mid<<=1){
            wn=complex(cos(pi/mid),type*sin(pi/mid));
            for(j=0;j<lim;j+=(mid<<1)){
                w=complex(1,0);
                for(k=0;k<mid;k++,w=w*wn){
                    x=a[j+k];y=a[j+k+mid]*w;
                    a[j+k]=x+y;
                    a[j+k+mid]=x-y;
                }
            }
        }
        if(type==-1) for(i=0;i<lim;i++) a[i].x=(ll)(a[i].x/lim+0.5);
    }
    void init(int len){
        memset(A,0,sizeof(A));memset(B,0,sizeof(B));
        cnt=0;lim=1;
        while(lim<=len) lim<<=1,cnt++;
        for(int i=0;i<lim;i++) r[i]=((r[i>>1]>>1)|((i&1)<<(cnt-1)));
    }
}
void mul(){
    using namespace FFT;
    fft(A,1);
    fft(B,1);
    for(int i=0;i<lim;i++) A[i]=A[i]*B[i];
    fft(A,-1);
}
double sqr(double x){
    return x*x;
}
ll n,k,x[20010],y[20010],pre1[20010],pre2[20010];
int main(){
    int T=read();
    while(T--){
        n=read();k=read();
        ll t1=0,t2=0,i,j,ans=1e15,tl,tm;
        
        for(i=0;i<k;i++){
            x[i]=read();
            t1+=(x[i]*x[i]);
            t2+=2*x[i];
        }
        for(i=0;i<k;i++){
            y[i]=read();
            y[i+k]=y[i]+n;
        }
        pre1[0]=y[0]*y[0];
        pre2[0]=y[0];
        for(i=1;i<(k<<1);i++){
            pre1[i]=pre1[i-1]+y[i]*y[i];
            pre2[i]=pre2[i-1]+y[i];
        }
        
        FFT::init(k*3);
        for(i=0;i<k;i++) FFT::A[i].x=x[k-1-i];
        for(i=0;i<(k<<1);i++) FFT::B[i].x=y[i];
        mul();
        
        for(i=k-1;i<(k<<1)-1;i++){
            j=i-k+1;
            tl=t1+pre1[j+k-1]-pre1[j-1];
            tm=t2-2*(pre2[j+k-1]-pre2[j-1]);//这里处理的是y的前缀和,以及y的平方的前缀和
            ll tmpc=(ll)(-(double)tm/(2.0*(double)k));

            ans=min(ans,(ll)(k*sqr(tmpc)+tmpc*tm+tl-2*FFT::A[i].x));//这里需要上下都摸♂索一下
            tmpc--;
            ans=min(ans,(ll)(k*sqr(tmpc)+tmpc*tm+tl-2*FFT::A[i].x));
            tmpc+=2;
            ans=min(ans,(ll)(k*sqr(tmpc)+tmpc*tm+tl-2*FFT::A[i].x));
        }
        printf("%lld\n",ans);
    }
}

转载于:https://www.cnblogs.com/dedicatus545/p/9722705.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值