【HNOI2017】礼物（FFT）

显然，\(y_i\) 加上 \(c\) 可以看成是 \(x_i\) 减去 \(c\)。

所以就变成了 \(x_i\) 加上一个整数（可正可负）。

现将 \(x\) 环拆成一个长度为 \(2n\) 的序列 \(a\)（复制一遍），把 \(y\) 环拆成一个长度为 \(n\) 的序列 \(b\)。

那么旋转操作就可以看成是 \(b\) 序列与 \(a\) 序列中每一个长度为 \(n\) 的子串匹配求值。

也就是说，求这个东西的最小值：\(\sum_{j=0}^{n-1}(a_{i+j}-b_j+c)^2\)（\(0\leq i< n\)）。

接下来推式子：（设 \(x\) 环的和是 \(suma\)，\(y\) 环的和是 \(sumb\)，\(x\) 环的平方和是 \(powa\)，\(y\) 环的平方和是 \(powa\)）

\[\begin{aligned} &\min_{i=0}^{n-1} \sum_{j=0}^{n-1}(a_{i+j}-b_j+c)^2\\ =&\min_{i=0}^{n-1}\sum_{j=0}^{n-1}(a_{i+j}-b_j)^2+c^2+2(a_{i+j}-b_j)c\\ =&\min_{i=0}^{n-1} [nc^2+\sum_{j=0}^{n-1}(a_{i+j}-b_j)^2+2c\sum_{j=0}^{n-1}(a_{i+j}-b_j)]\\ =&\min_{i=0}^{n-1} [nc^2+\sum_{j=0}^{n-1}(a_{i+j}^2+b_j^2-2a_{i+j}b_j)+2c(\sum_{j=0}^{n-1} a_{i+j}-\sum_{j=0}^{n-1} b_j)]\\ =&\min_{i=0}^{n-1} (nc^2+powa+powb-2\sum_{j=0}^{n-1}a_{i+j}b_j+2c(suma-sumb)) \end{aligned}\]

发现 \(powa+powb\) 是定值，而和 \(c\) 有关的 \(nc^2+2c(suma-sumb)\) 可以用二次函数最值 \(O(1)\) 求，也就是说我们只需要求 \(-2\sum_{j=0}^{n-1}a_{i+j}b_j\) 的最小值，即 \(\sum_{j=0}^{n-1}a_{i+j}b_j\) 的最大值。

我们可以把这个形式改变一下，把 \(\sum_{j=0}^{n-1}a_{i+j}b_j\) 改成 \(\sum_{i-j=k}a_ib_j\)（\(k\) 是旋转角度）。

令 \(S(k)=\sum_{i-j=k}a_ib_j\)，设 \(\widehat{a_i}=a_{2n-i}\)，则：

\[S(k)=\sum_{i-j=k}^{0\leq j< n}a_ib_j=\sum_{(2n-i)-j=k}^{0\leq j< n}\widehat{a_i}b_j=\sum_{i+j=2n-k}^{0\leq j< n}\widehat{a_i}b_j \]

把 \(n-k\) 代入：

\[S(2n-k)=\sum_{i+j=2n-(2n-k)}^{0\leq j< n}\widehat{a_i}b_j=\sum_{i+j=k}^{0\leq j< n}\widehat{a_i}b_j \]

有木有觉得很熟悉？

我们可以把 \(\widehat{a_i}\) 看成一个多项式 \(A\) 的系数，\(b_j\) 看成另一个多项式 \(B\) 的系数，那么 \(S(2n-k)\) 就是 \(A\times B\) 第 \(k\) 项的系数。

看会我们原来的问题：求 \(\sum_{i-j=k}a_ib_j\) 的最大值，也就是求 \(S(k)\) 的最大值（\(0\leq k<n\)），也就是 \(A\times B\) 第 \((n+1)\sim 2n\) 项系数的最大值。

那么现在就好处理了，我们先用 FFT 把 \(A\times B\) 算出来，再取最大值。

代码如下：

#include<bits/stdc++.h>

#define N 50010
#define PN 262200
#define INF 0x7fffffff

using namespace std;

struct Complex
{
	double x,y;
	Complex(){};
	Complex(double xx,double yy){x=xx,y=yy;}
}a[PN],b[PN],c[PN];

Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

const double pi=acos(-1);

int n,m,suma,sumb,powa,powb;
int limit=1,bit,rev[PN];

void FFT(Complex *a,int opt)
{
	for(int i=0;i<limit;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		Complex wn=Complex(cos(pi/mid),opt*sin(pi/mid));
		for(int i=0,len=(mid<<1);i<limit;i+=len)
		{
			Complex now=Complex(1,0);
			for(int j=0;j<mid;j++,now=now*wn)
			{
				Complex x=a[i+j],y=now*a[i+mid+j];
				a[i+j]=x+y,a[i+mid+j]=x-y;
			}
		}
	}
	if(opt==-1)
		for(int i=0;i<limit;i++)
			a[i].x/=limit;
}

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<n;i++)
	{
		scanf("%lf",&a[i].x);
		a[n+i].x=a[i].x;
		suma+=a[i].x;
		powa+=a[i].x*a[i].x;
	}
	for(int i=0;i<n;i++)
	{
		scanf("%lf",&b[i].x);
		sumb+=b[i].x;
		powb+=b[i].x*b[i].x;
	}
	int d=round(-1.0*(suma-sumb)/n);
	reverse(a,a+2*n+1);//求a^
	while(limit<=3*n)
		limit<<=1,bit++;
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	FFT(a,1),FFT(b,1);
	for(int i=0;i<limit;i++)
		c[i]=a[i]*b[i];
	FFT(c,-1);
	int ans=-INF;
	for(int i=n+1;i<=2*n;i++)
		ans=max(ans,(int)round(c[i].x));//这里用round是怕被卡精度
	printf("%d\n",n*d*d+2*d*(suma-sumb)+powa+powb-2*ans);
	return 0;
}