// author:latstars
// time:26/092016 19:52
// Problem:hihoCoder #1388 : Periodic Signal
// Class:FFT plus special judge
// or NTT + CRT
// 题意:求最小的 sum(A[i]-B[i+k% n])^2 for i=0 to n-1 for k=0 to n-1
// 思路:展开 得到sum (A[i]^2+B[i+k %n]^2-2*A[i]*B[i+k %n]) for i=0 to n-1
// 前面的A[i]^2 和 B[i+k %n]^2 的sum 是一个定值
// 那么只要求sum(A[i]*B[i+k%n]) for i=0 to n-1的最大值就可以
// 是一个卷积
// 构造多项式乘法
// A[0] A[1] A[2] ...... A[n-1] 0 0 0 ...... 0 0 0
// B[n-1] B[n-2] ...... B[0] B[n-1] B[n-2] ...... B[2] B[1] B[0]
// 所以 呢 C[n-1] =A[0]B[0]+A[1]B[1] ....
// C[n]=A[0]B[n-1]+A[1]B[0]+...
// C[2*n-2]=A[0]B[1]+A[1]B[2]+...
// 但是会出现精度误差,那么可以这样找出最大的卷积对应的位置
// 然后再算一下原序列
// 另外一种方法是用NTT 但是不能使用p为long long 的模数,因为会暴long long
// 所以可以选择多个不同的模数,然后使用中国剩余定理还原原来的数
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<complex>
#include<vector>
#include<algorithm>
using namespace std;
const int maxn=(1<<21)+1000;
const int maxm=(1<<18)+10;
const double PI=acos(-1.0);
typedef complex<long double> Complex;
void FFT(Complex * a,int n,int inv)
{
for(int i=0,j=0;i<n;i++){
if(j>i) swap(a[i],a[j]);
int k=n;
while(j&(k>>=1)) j&=~k;
j|=k;
}
for(int s=1;s<n;s<<=1){//s=2^0 ,s wei chuli de DFT yuanshu de yiban
double ap=PI/s*inv;//ap=2*PI/(s*2)*inv
Complex wn=Complex(1.0,0.0);
Complex w=Complex(cos(ap),sin(ap));//w=
for(int k=0;k<s;k++){
// cout<<wn<<endl;
// w=Complex(cos(ap*k),sin(ap*k));
for(int e=k;e<n;e+=s<<1){
int o=e+s;
Complex u=wn*a[o];
a[o]=a[e]-u;
a[e]+=u;
}
wn=wn*w;
}
}
if(inv==-1){
for(int i=0;i<n;i++) a[i]/=n;
}
}
typedef long long LL;
LL A[maxm],B[maxm];
Complex C[maxn];
Complex D[maxn];
int T;
int n;
void solve(void)
{
}
int main()
{
scanf("%d",&T);
while(T--){
LL ans=0;
memset(A,0,sizeof(A));
memset(B,0,sizeof(B));
scanf("%d",&n);
for(int i=0;i<n;i++){
scanf("%lld",&A[i]);
// A[i]=1<<20;
ans=ans+A[i]*A[i];
// cout<<ans<<endl;
}
for(int i=0;i<n;i++){
scanf("%lld",&B[n-i-1]);
// B[n-i-1]=1<<20;
ans=ans+B[n-i-1]*B[n-i-1];
}
for(int i=n;i<2*n;i++){
A[i]=0;
B[i]=B[i-n];
}
n=n*2;
int len=1;
while(len<=n) len<<=1;
len<<=1;
// cout<<len<<endl;
// cout<<ans<<endl;
for(int i=0;i<n;i++){
C[i]=Complex(A[i],0.0);
D[i]=Complex(B[i],0.0);
}
for(int i=n;i<len;i++){
C[i]=D[i]=Complex(0.0,0.0);
}
FFT(C,len,1);
FFT(D,len,1);
for(int i=0;i<len;i++){
C[i]=C[i]*D[i];
}
FFT(C,len,-1);
LL tmp=0;
int pos=0;
for(int i=n/2;i<n;i++){//
if(tmp<=((C[i].real()+0.1))){
tmp=(C[i].real()+0.1);
pos=i;
}
// cout<<i<<" "<<(LL)(C[i].real()+0.1)<<endl;
}
tmp=0;
for(int i=0;i<n/2;i++){
tmp+=A[i]*B[pos-i];
}
printf("%lld\n",ans-2*tmp);
}
return 0;
}
2016ICPC北京站网络预选赛F hihoCoder1388 Periodic Signal,FFT|NTT+CRT
最新推荐文章于 2024-09-23 09:35:15 发布