【hihocoder 1388】 【NTT或者FFT 循环矩阵】

传送门：点击打开链接

题意：

给你两个序列，求一个差值平方和最小
具体求法就是
min{(a1-b1)²+...+(an-bn)²,(a1-b2)²+...+(an-b1)²，...，(a1-bn)²+...+(an-b1)²}

思路：

化简的得
∑ni=1a2i+∑ni=1b2i−2(ai∗bj)(iϵn,jϵn)
现在问题就变成了怎么快速求出ai∗bj 的最大值
普通方法的时间复杂度O(n2),肯定要超时。
考虑下面的循环矩阵

⎡⎣⎢⎢⎢⎢⎢⎢a1anan−1.a2a2a1an.a3a3a2a1.a4...............anan−1an−2.a1⎤⎦⎥⎥⎥⎥⎥⎥∗⎡⎣⎢⎢⎢⎢⎢⎢b1b2b3.bn⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢c1c2c3.cn⎤⎦⎥⎥⎥⎥⎥⎥

只需要求cn 的最大值即可。
可以参考下面一个链接
http://www.zybang.com/question/dd8b336ad690b3e2f9aa4dc0b596e1ea.html
将循环矩阵*向量转换为求卷积
所以用快速傅里叶变换将时间复杂度降为了O(nlog2n)
具体做法是
令a→=(a1,a2...,an)
令b→=(bn,bn−1,...b2,b1,bn,bn−1...b2)
之后对a→ 求 c→=DFT(a→)
对b→ 求d→= DFT(b→)
时间复杂度为O(nlog2n)
然后求S→=c→∗d→
时间复杂度为O(n)
再对S→ 求ans−→−=IDFT(S→)
时间复杂度为O(nlog2n)
最后求ans−→− 中最大的元素即可 

用FFT要注意精度问题，一个无赖的办法是根据FFT 的结果求 K，然后再自己算一遍得到最后答案，具体代码见Here

FFT文章推荐：http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15

和http://blog.csdn.net/acdreamers/article/details/39005227

NTT文章推荐：https://riteme.github.io/blog/2016-8-22/ntt.html

和http://m.blog.csdn.net/article/details?id=39026505

NNT中的费马素数表：http://blog.miskcoo.com/2014/07/fft-prime-table

另一种方法是用NTT（快速数论变换）算完直接比较出最大值，long long版的板子见Here

代码：

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;#define MAXN 280000 const long long P=50000000001507329LL; // 190734863287 * 2 ^ 18 + 1//const int P=1004535809; // 479 * 2 ^ 21 + 1//const int P=998244353; // 119 * 2 ^ 23 + 1const int G=3; long long mul(long long x,long long y){  return (x*y-(long long)(x/(long double)P*y+1e-3)*P+P)%P;}long long qpow(long long x,long long k,long long p){  long long ret=1;  while(k){    if(k&1) ret=mul(ret,x);    k>>=1;    x=mul(x,x);  }  return ret;} long long wn[25];void getwn(){  for(int i=1; i<=18; ++i){    int t=1<<i;    wn[i]=qpow(G,(P-1)/t,P);  }} int len;void NTT(long long y[],int op){  for(int i=1,j=len>>1,k; i<len-1; ++i){    if(i<j) swap(y[i],y[j]);    k=len>>1;    while(j>=k){      j-=k;      k>>=1;    }    if(j<k) j+=k;  }  int id=0;  for(int h=2; h<=len; h<<=1) {    ++id;    for(int i=0; i<len; i+=h){      long long w=1;      for(int j=i; j<i+(h>>1); ++j){        long long u=y[j],t=mul(y[j+h/2],w);        y[j]=u+t;        if(y[j]>=P) y[j]-=P;        y[j+h/2]=u-t+P;        if(y[j+h/2]>=P) y[j+h/2]-=P;        w=mul(w,wn[id]);      }    }  }  if(op==-1){  for(int i=1; i<len/2; ++i) swap(y[i],y[len-i]);  long long inv=qpow(len,P-2,P);  for(int i=0; i<len; ++i) y[i]=mul(y[i],inv);  }}void Convolution(long long A[],long long B[],int n){  for(len=1; len<(n<<1); len<<=1);  for(int i=n; i<len; ++i){    A[i]=B[i]=0;  }    NTT(A,1); NTT(B,1);  for(int i=0; i<len; ++i){    A[i]=mul(A[i],B[i]);  }  NTT(A,-1);} long long A[MAXN],B[MAXN];int main(){  getwn();  int t,n;  scanf("%d",&t);  while(t--){    scanf("%d",&n);    long long ans=0;    for(int i=0; i<n; ++i){      scanf("%lld",&A[i]);      ans+=A[i]*A[i];    }    for(int i=0; i<n; ++i){      scanf("%lld",&B[n-i-1]);      ans+=B[n-i-1]*B[n-i-1];    }    for(int i=0; i<n; ++i){      A[i+n]=A[i];      B[i+n]=0;    }    Convolution(A,B,2*n);    long long mx=0;    for(int i=n; i<2*n; ++i){      mx=max(mx,A[i]);    }    printf("%lld\n",ans-2*mx);  }  return 0;}

常用费马素数和原根的表：

/*是这样的，这几天在写 FFT，由于是在模意义下的，需要各种素数……然后就打了个表方便以后查了如果 r*2^k+1 是个素数，那么在mod r*2^k+1意义下，可以处理 2^k 以内规模的数据，2281701377=17*2^27+1  是一个挺好的数，平方刚好不会爆 long long1004535809=479*2^21+1 加起来刚好不会爆 int 也不错下面是刚刚打出来的表格（gg 是mod(r*2^k+1)的原根）素数  rr  kk  gg3   1   1   25   1   2   217  1   4   397  3   5   5193 3   6   5257 1   8   37681    15  9   1712289   3   12  1140961   5   13  365537   1   16  3786433  3   18  105767169 11  19  37340033 7   20  323068673    11  21  3104857601   25  22  3167772161   5   25  3469762049   7   26  31004535809  479 21  32013265921  15  27  312281701377  17  27  33221225473  3   30  575161927681 35  31  377309411329 9   33  7206158430209    3   36  222061584302081   15  37  72748779069441   5   39  36597069766657   3   41  539582418599937  9   42  579164837199873  9   43  5263882790666241 15  44  71231453023109121    35  45  31337006139375617    19  46  33799912185593857    27  47  54222124650659841    15  48  197881299347898369    7   50  631525197391593473   7   52  3180143985094819841  5   55  61945555039024054273 27  56  54179340454199820289 29  57  3*/