【hihoCoder #1388】【2016 acm 北京网络赛 F题】【FFT】 Periodic Signal

传送门：hihocoder 1388

时间限制: 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

思路：

scrolling="no" style="border-width: medium; border-style: none; width: 686px; height: 109px;" src="data:text/html;charset=utf8,%3Cimg%20id=%22img%22%20src=%22https://media.hihocoder.com/problem_images/20160923icpc/formulauysyuhjs023.png?_=0.5634232694209171%22%20style=%22border:none;max-width:686px%22%3E%3Cscript%3Ewindow.onload%20=%20function%20%28%29%20%7Bvar%20img%20=%20document.getElementById%28%27img%27%29;%20window.parent.postMessage%28%7BiframeId:%27iframe_0.9635411305986604%27,width:img.width,height:img.height%7D,%20%27http://www.cnblogs.com%27%29;%7D%3C/script%3E" id="iframe_0.9635411305986604" frameborder="0">

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

 

即求A[]与B[]的循环卷积。 FFT求解。

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

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <math.h>
#define pl(x) cout << #x << "= " << x << endl;
using namespace std;
typedef long long ll;

const double PI = acos(-1.0);
//复数结构体
struct complex{
  double r,i;
  complex(double _r = 0.0,double _i = 0.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);
  }
};
/*
 * 进行FFT和IFFT前的反转变换。
 * 位置i和 （i二进制反转后位置）互换
 * len必须去2的幂
 */
void change(complex y[],int len){
  int i,j,k;
  for(i = 1, j = len/2;i < len-1; i++){
    if(i < j)swap(y[i],y[j]);
    //交换互为小标反转的元素，i<j保证交换一次
    //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
    k = len/2;
    while( j >= k){
      j -= k;
      k /= 2;
    }
     if(j < k) j += k;
  }
}
/*
 * 做FFT
 * len必须为2^k形式，
 * on==1时是DFT，on==-1时是IDFT
 */
void fft(complex y[],int len,int on){
  change(y,len);
  for(int h = 2; h <= len; h <<= 1){
    complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
    for(int j = 0;j < len;j+=h){
      complex w(1,0);
      for(int k = j;k < j+h/2;k++){
        complex u = y[k];
        complex t = w*y[k+h/2];
        y[k] = u+t;
        y[k+h/2] = u-t;
        w = w*wn;
      }
    }
  }
  if(on == -1)
  for(int i = 0;i < len;i++)
    y[i].r /= len;
}

const int maxn = 240040;
complex x1[maxn],x2[maxn];
ll a[maxn/4],b[maxn/4];  //原数组
ll sum[maxn];  //FFT结果

void init(){
	memset(sum, 0, sizeof(sum));
	memset(x1, 0, sizeof(x1));
	memset(x2, 0, sizeof(x2));
}

int main(){
	int T;
	scanf("%d",&T);
  while(T--){
/*    int len1 = strlen(str1);
    int len2 = strlen(str2);*/
    int n;ll sa=0,sb=0;
    init();
    scanf("%d",&n);
    for(int i=0; i<n; i++){ scanf("%lld",&a[i]); sa+=a[i]*a[i]; }
    for(int i=0; i<n; i++){ scanf("%lld",&b[i]); sb+=b[i]*b[i]; }

    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 = len1;i < len;i++)
      x1[i] = complex(0,0);*/
    for(int i=0; i<n; i++)
      x2[i] = complex(b[n-i-1], 0);//反向
/*    for(int i = len2;i < len;i++)
      x2[i] = complex(0,0);*/
    //求DFT
    fft(x1,len,1);
    fft(x2,len,1);
    for(int i = 0;i < len;i++)
      x1[i] = x1[i]*x2[i];
    fft(x1,len,-1);
    for(int i = 0;i < len;i++)
      sum[i] = (ll)(x1[i].r+0.5);
    ll res=sum[n-1]+sum[2*n-1];int k=0;
    for(int i=0; i<=n-2; i++){
    	if(res<sum[i]+sum[i+n]){ //可以写两个四项式相乘的例子，sum[i]+sum[i+n]就是循环卷积和
    		res=sum[i]+sum[i+n];//注意此时得到的res会有很小的浮点精度误差
    		k=n-i-1;//但是记录的k,这个是正确的
    	}
    }	
    res=0;
    for(int i=0; i<n; i++){
    	res+=a[i]*b[(i+k)%n]; //重新算一遍得到最后答案
    }
    ll ans=sa+sb-2*res;
    printf("%lld\n",ans);
  }  
  return 0;	
}