hihocoder 1388 2016北京网络赛1006 Periodic Signal（卡精度的FFT）

#1388 : Periodic Signal

时间限制:5000ms

单点时限:5000ms

内存限制:256MB

描述

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.

样例输入

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

题意：给定一个  A  数组和  B  数组，求给定的这个函数值的最小值。

题解：

A[]的平方和 与 B[]的平方和可以直接求出。所以只要求出的最大值即可得到答案。

很明显是求求A[]与B[]的卷积。 FFT搞定。

注意由于数据较大，FFT会出现精度问题。最后结果会有浮点精度误差，但是由结果得到的 k 是正确的，所以唯一的办法是根据FFT的结果求 K，然后再自己算一遍得到最后答案。

注：题解的标准做法是找两个 10910^910​9​​ 左右模数 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;
}