hihocoder 1388 Periodic Signal

Periodic Signal

时间限制:5000ms

单点时限:5000ms

内存限制:256MB

Description

Profess X is an expert in signal processing. He has a device which can send a particular 1 second signal repeatedly. The signal is A0 ... An-1 under n Hz sampling.

One day, the device fell on the ground accidentally. Profess X wanted to check whether the device can still work properly. So he ran another n Hz sampling to the fallen device and got B0 ... Bn-1.

To compare two periodic signals, Profess X define the DIFFERENCE of signal A and B as follow:

You may assume that two signals are the same if their DIFFERENCE is small enough.
Profess X is too busy to calculate this value. So the calculation is on you.

Input

The first line contains a single integer T, indicating the number of test cases.

In each test case, the first line contains an integer n. The second line contains n integers, A0 ... An-1. The third line contains n integers, B0 ... Bn-1.

T≤40 including several small test cases and no more than 4 large test cases.

For small test cases, 0<n≤6⋅103.

For large test cases, 0<n≤6⋅104.

For all test cases, 0≤Ai,Bi<220.

Output

For each test case, print the answer in a single line.

Sample Input

2
9
3 0 1 4 1 5 9 2 6
5 3 5 8 9 7 9 3 2
5
1 2 3 4 5
2 3 4 5 1

Sample Output

80
0

/** hihocoder 1388 Periodic Signal** 把式子变形一下就是求Ai*Bi+k求和的最大值，想到用FFT来求* 会因为精度问题不能过，这是找出最小的那个k，然后重新算即可。*/#include <bits/stdc++.h>
using namespace std;const int MAXN = 200000+10;
const double PI = acos(-1.0);struct Complex
{double r,i;Complex(double _r=0,double _i=0):r(_r),i(_i){}Complex operator + (const Complex& rhs){return Complex(r+rhs.r,i+rhs.i);}Complex operator - (const Complex& rhs){return Complex(r-rhs.r,i-rhs.i);}Complex operator * (const Complex &rhs){return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i);}
};
/** 进行FFT和IFFT前的反转变换。* 位置i和 （i二进制反转后位置）互换* len必须取2的幂*/
void Rader(Complex F[],int len)
{int j = len >> 1;for(int i = 1;i < len - 1;++i){if(i < j) swap(F[i],F[j]);  // reverseint k = len>>1;while(j>=k){j -= k;k >>= 1;}if(j < k) j += k;}
}
/** 做FFT* len必须为2^k形式，* on==1时是DFT，on==-1时是IDFT*/
void FFT(Complex F[],int len,int t)
{Rader(F,len);for(int h=2;h<=len;h<<=1){Complex wn(cos(-t*2*PI/h),sin(-t*2*PI/h));for(int j=0;j<len;j+=h){Complex E(1,0); //旋转因子for(int k=j;k<j+h/2;++k){Complex u = F[k];Complex v = E*F[k+h/2];F[k] = u+v;F[k+h/2] = u-v;E=E*wn;}}}if(t==-1)   //IDFTfor(int i=0;i<len;++i)F[i].r/=len;
}
void Conv(Complex a[],Complex b[],int len) //求卷积
{FFT(a,len,1);FFT(b,len,1);for(int i=0;i<len;++i) a[i] = a[i]*b[i];FFT(a,len,-1);
}
Complex va[MAXN],vb[MAXN];
int A[MAXN],B[MAXN],n,len;
int AA[MAXN];
void init()
{scanf("%d",&n);for(int i=0;i<n;i++) scanf("%d",&A[i]);for(int i=0;i<n;i++) scanf("%d",&B[i]);for(int i=n;i<2*n;i++) B[i]=B[i-n];for(int i=0;i<n;i++) AA[i]=A[n-i-1];len=1;while(len<2*n) len<<=1;for(int i=0;i<n;i++) va[i]=Complex(AA[i],0);for(int i=n;i<len;i++) va[i]=Complex(0,0);for(int i=0;i<2*n;i++) vb[i]=Complex(B[i],0);for(int i=2*n;i<len;i++) vb[i]=Complex(0,0);
}
void solve()
{init();Conv(va,vb,len);double ans=0;int k;for(int i=n-1;i<=2*n-2;i++){if(va[i].r>ans){ans=va[i].r;k=i;}}k-=n-1;long long a=0,b=0,ab=0;for(int i=0;i<n;i++) a+=(long long)A[i]*A[i];for(int i=0;i<n;i++) b+=(long long)B[i]*B[i];for(int i=0;i<n;i++) ab+=(long long)A[i]*B[i+k];long long res=a+b-2*ab;printf("%lld\n",res);
}
int main()
{int T;int n;scanf("%d",&T);while(T--){solve();}return 0;
}

转载于:https://www.cnblogs.com/wangdongkai/p/5906954.html