【BZOJ4827】礼物（AHOI&HNOI2017）-FFT

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

#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;
}