#1388 : Periodic Signal
-
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
样例输出
-
80 0
描述
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.
输入
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.
输出
For each test case, print the answer in a single line.
题意:给定一个 A 数组和 B 数组,求给定的这个函数值的最小值。
题解:
A[]的平方和 与 B[]的平方和可以直接求出。所以只要求出的最大值即可得到答案。
很明显是求求A[]与B[]的卷积。 FFT搞定。
注意由于数据较大,FFT会出现精度问题。最后结果会有浮点精度误差,但是由结果得到的 k 是正确的,所以唯一的办法是根据FFT的结果求 K,然后再自己算一遍得到最后答案。
注:题解的标准做法是找两个 10910^9109
左右模数 NTT 后 CRT 。。。。
AC代码:
#include<bits/stdc++.h>
using namespace std;
const long double pi = acos(-1.0);
typedef long long ll;
struct Complex // 复数类结构
{
long double r,i;
Complex(long double _r = 0,long double _i = 0)
{
r = _r; i = _i;
}
Complex operator +(const Complex &b)
{
return Complex(r+b.r,i+b.i);
}
Complex operator -(const Complex &b)
{
return Complex(r-b.r,i-b.i);
}
Complex operator *(const Complex &b)
{
return Complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};
//雷德算法--倒位序
void rader(Complex F[],int len) 二进制平摊反转置换:O(logn) :len = 2^M,reverse F[i] with F[j] j为i二进制反转
{
int j = len >> 1;
for(int i = 1;i < len - 1;++i)
{
if(i < j) swap(F[i],F[j]); // reverse
int k = len>>1;
while(j>=k)
{
j -= k;
k >>= 1;
}
if(j < k) j += k;
}
}
//FFT实现
void FFT(Complex F[],int len,int t)
{
rader(F,len);
for(int h=2;h<=len;h<<=1) //分治后计算长度为h的DFT
{
Complex wn(cos(-t*2*pi/h),sin(-t*2*pi/h)); //单位复根e^(2*PI/m)用欧拉公式展开
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) //IDFT
for(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);
}
const int MAXN = 240040;
Complex x1[MAXN],x2[MAXN];
ll a[MAXN/4],b[MAXN/4]; //原数组
ll num[MAXN]; //FFT结果
int n;
ll suma,sumb;
void solve()//把A加长一倍,B反转,构造多项式,这样多项式乘积指数[n,2n)的系数就是各个位置的结果了
{
int len = 1;
while( len < 2*n ) len <<= 1;
for(int i = 0;i < n;i++)
{
x1[i] = Complex(a[i],0);
}
for(int i=n;i<len;i++)
{
x1[i]=Complex(0,0);
}
for(int i = 0;i < n;i++)
{
x2[i] = Complex(b[n-i-1],0);
}
Conv(x1,x2,len);
for(int i = 0;i < len;i++)
{
num[i] = (ll)(x1[i].r+0.5);
}
}
void gao()
{
ll ret=num[n-1];
int k=0;
for(int i=0;i<n-2;i++)
{
if(ret<num[i]+num[i+n])
{
ret=num[i]+num[i+n];
k=n-1-i;
}
//注意:此时得到的ret会有很小的浮点精度误差,
}
ret=0;
for(int i=0;i<n;i++)
{
ret+=a[i]*b[(i+k)%n]; //重新算一遍得到最后答案
}
ll ans=suma+sumb-2*ret;
printf("%lld\n",ans);
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
suma=0;
sumb=0;
memset(num,0,sizeof(num));
memset(x1,0,sizeof(x1));
memset(x2,0,sizeof(x2));
scanf("%d",&n);
for(int i = 0;i < n;i++)
{
scanf("%lld",&a[i]);
suma+=a[i]*a[i];
}
for(int i = 0;i < n;i++)
{
scanf("%lld",&b[i]);
sumb+=b[i]*b[i];
}
solve();
gao();
}
return 0;
}