【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 ),肯定要超时。
考虑下面的循环矩阵

只需要求 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 + 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("%lldn",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
*/

最后

以上就是暴躁火龙果为你收集整理的【hihocoder 1388】 【NTT或者FFT 循环矩阵】的全部内容，希望文章能够帮你解决【hihocoder 1388】 【NTT或者FFT 循环矩阵】所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。