显然,\(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;
}
标签:sum,FFT,widehat,leq,Complex,HNOI2017,2n,PN,礼物
From: https://www.cnblogs.com/ez-lcw/p/16837452.html