HihoCoder 1388 Periodic Signal 循环卷积(FFT)

题意:给出n个数的序列a,b ,n<=6e4.

拆来后得: tot - 2*max( a[i]*b[(i+k)%n] i=0~n-1 )k=0,..n-1)
先写成矩阵形式后,也就是要求这n个乘积累加和中的最大值.
循环卷积的构造方式:点击打开链接
构造两个多项式,a[0]+a[1]*x+a[2]*x^2+..a[n-1]*x^(n-1).
b[n-1]+b[n-2]*x+.....b[0]*x^(n-1) + b[n-1]*x^n+....b[0]*x^(2n-1)
将b逆序后重新写两遍,和a相乘 则相乘后X^(n-1)的系数就是k==0时结果.其余类似.
fft爆精度 求最大系数位置 然后代入答案求一下即可.

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 150005
#define pi acos(-1.0)
using namespace std;
int l,T,n,num,wz;
long long ans,tot,sum[N],A[N],B[N];
struct complex
{  
    double r,i;  
    complex(double real=0.0,double image=0.0){r=real;i=image;}
    complex operator + (const complex o)
    {
        return complex(r+o.r,i+o.i);
    }
    complex operator - (const complex o)
    {
        return complex(r-o.r,i-o.i);
    }
    complex operator * (const complex o)
    {
        return complex(r*o.r-i*o.i,r*o.i+i*o.r);
    }  
}x1[N],x2[N];
void brc(complex *y,int l)
{  
    register int i,j,k;
    for(i=1,j=l/2;i<l-1;i++)  
    {
        if(i<j)swap(y[i],y[j]);
        k=l/2;
        while(j>=k)j-=k,k/=2;
        if(j<k)j+=k;
    }
}
void fft(complex *y,int l,double on)
{
    register int h,i,j,k;
    complex u,t;
    brc(y,l);
    for(h=2;h<=l;h<<=1)
    {  
        complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(j=0;j<l;j+=h)
        {
            complex w(1,0);
            for(k=j;k<j+h/2;k++)
            {
                u=y[k];
                t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)for(i=0;i<l;i++)y[i].r/=l;
}
long long Getans(int k)
{
    ans=0;
    for(int i=0;i<n;i++)ans+=(A[i]-B[(i+k)%n])*(A[i]-B[(i+k)%n]);
    return ans;
}
int main()  
{
    register int i,j;
    scanf("%d",&T);
    while(T--)
    {
    	l=1;
    	tot=0;
    	scanf("%d",&n);
    	while(l<2*n) l<<=1;
    	for(int i=n-1;i>=0;i--)
    	{
    		scanf("%d",&num);
    		x1[i].r=num,x1[i].i=0;
    		A[n-i-1]=num;
		}
		for(int i=n;i<2*n-1;i++)
			x1[i]=x1[i-n];
		for(int i=2*n-1;i<l;i++)
			x1[i].r=x1[i].i=0;
		
		for(int i=0;i<n;i++)
		{
			scanf("%d",&num);
			x2[i].r=num,x2[i].i=0;
			B[i]=num;
		}
		for(int i=n;i<l;i++)
			x2[i].r=x2[i].i=0;
		fft(x1,l,1),fft(x2,l,1);//DFT.
		
		for(int i=0;i<l;i++)//mul
			x1[i]=x1[i]*x2[i];
			
		fft(x1,l,-1);//IDFT
		
		double m=-1;
		int pos;
		for(int i=n-1;i<2*n;i++)
			if(m<x1[i].r)
				m=x1[i].r,pos=i;
		printf("%lld\n",Getans((pos+1-n)));
	}
    return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值