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

最新推荐文章于 2024-11-13 23:40:34 发布

原创最新推荐文章于 2024-11-13 23:40:34 发布 · 1.1k 阅读

0 ·

CC 4.0 BY-SA版权

online judge hihocoder 同时被 2 个专栏收录

3 篇文章

订阅专栏

------FFT&&NTT

3 篇文章

订阅专栏

本文介绍如何使用快速傅立叶变换（FFT）来高效解决两序列间最小差值平方和的问题，并提供了一种利用循环矩阵和卷积的方法，通过C++代码实现了这一解决方案。

传送门：点击打开链接

题意：

给你两个序列，求一个差值平方和最小
具体求法就是
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 ),肯定要超时。
考虑下面的循环矩阵

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ a 1 a n a n - 1 . a 2 a 2 a 1 a n . a 3 a 3 a 2 a 1 . a 4 . . . . . . . . . . . . . . . a n a n - 1 a n - 2 . a 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ * ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ b 1 b 2 b 3 . b n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ c 1 c 2 c 3 . c n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

只需要求

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

和http://blog.youkuaiyun.com/acdreamers/article/details/39005227

和http://m.blog.youkuaiyun.com/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 + 1
const 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 long

1004535809=479*2^21+1 加起来刚好不会爆 int 也不错

下面是刚刚打出来的表格（gg 是mod(r*2^k+1)的原根）


素数  rr  kk  gg
3   1   1   2
5   1   2   2
17  1   4   3
97  3   5   5
193 3   6   5
257 1   8   3
7681    15  9   17
12289   3   12  11
40961   5   13  3
65537   1   16  3
786433  3   18  10
5767169 11  19  3
7340033 7   20  3
23068673    11  21  3
104857601   25  22  3
167772161   5   25  3
469762049   7   26  3
1004535809  479 21  3
2013265921  15  27  31
2281701377  17  27  3
3221225473  3   30  5
75161927681 35  31  3
77309411329 9   33  7
206158430209    3   36  22
2061584302081   15  37  7
2748779069441   5   39  3
6597069766657   3   41  5
39582418599937  9   42  5
79164837199873  9   43  5
263882790666241 15  44  7
1231453023109121    35  45  3
1337006139375617    19  46  3
3799912185593857    27  47  5
4222124650659841    15  48  19
7881299347898369    7   50  6
31525197391593473   7   52  3
180143985094819841  5   55  6
1945555039024054273 27  56  5
4179340454199820289 29  57  3
*/