hihoCoder1388 Periodic Signal（2016北京网赛F：NTT）

题目

Source

http://hihocoder.com/problemset/problem/1388

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*10³.

For large test cases, 0<n≤6*10⁴.

For all test cases, 0≤Ai,Bi<2²⁰.

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

分析

题目大概说k多少时上面的式子值最小。

展开那个式子可以得到，$sum A_i^2 + sum B_i^2 - 2sum A_iB_{i+k}$。前面是固定的，问题相当于最小化后面那一部分。

容易发现这是个可以用FFT快速求出的卷积形式。。

把A加长一倍，B反转，构造多项式，这样多项式乘积指数[n,2n)的系数就是各个位置的结果了。

直接FFT是不行的，千百亿的浮点数运算精度不够，比赛时搞了30多发没搞出来。。

我用了个NTT的板子A了，http://blog.csdn.net/u013654696/article/details/48653057，里面除了附赠了几个模数外还有一个O(1)快速乘不明觉厉。

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define MAXN 262144

const long long P=50000000001507329LL; // 190734863287 * 2 ^ 18 + 1
//const int P=1004535809; // 479 * 2 ^ 21 + 1
//const int P=998244353; // 119 * 2 ^ 23 + 1
const int G=3;

long long mul(long long x,long long y){
	return (x*y-(long long)(x/(long double)P*y+1e-3)*P+P)%P;
}
long long qpow(long long x,long long k,long long p){
	long long ret=1;
	while(k){
		if(k&1) ret=mul(ret,x);
		k>>=1;
		x=mul(x,x);
	}
	return ret;
}

long long wn[25];
void getwn(){
	for(int i=1; i<=18; ++i){
		int t=1<<i;
		wn[i]=qpow(G,(P-1)/t,P);
	}
}

int len;
void NTT(long long y[],int op){
	for(int i=1,j=len>>1,k; i<len-1; ++i){
		if(i<j) swap(y[i],y[j]);
		k=len>>1;
		while(j>=k){
			j-=k;
			k>>=1;
		}
		if(j<k) j+=k;
	}
	int id=0;
	for(int h=2; h<=len; h<<=1) {
		++id;
		for(int i=0; i<len; i+=h){
			long long w=1;
			for(int j=i; j<i+(h>>1); ++j){
				long long u=y[j],t=mul(y[j+h/2],w);
				y[j]=u+t;
				if(y[j]>=P) y[j]-=P;
				y[j+h/2]=u-t+P;
				if(y[j+h/2]>=P) y[j+h/2]-=P;
				w=mul(w,wn[id]);
			}
		}
    }
    if(op==-1){
		for(int i=1; i<len/2; ++i) swap(y[i],y[len-i]);
		long long inv=qpow(len,P-2,P);
		for(int i=0; i<len; ++i) y[i]=mul(y[i],inv);
    }
}
void Convolution(long long A[],long long B[],int n){
	for(len=1; len<(n<<1); len<<=1);
	for(int i=n; i<len; ++i){
		A[i]=B[i]=0;
	}
	
	NTT(A,1); NTT(B,1);
	for(int i=0; i<len; ++i){
		A[i]=mul(A[i],B[i]);
	}
	NTT(A,-1);
}

long long A[MAXN],B[MAXN];
int main(){
	getwn();
	int t,n;
	scanf("%d",&t);
	while(t--){
		scanf("%d",&n);
		long long ans=0;
		for(int i=0; i<n; ++i){
			scanf("%lld",&A[i]);
			ans+=A[i]*A[i];
		}
		for(int i=0; i<n; ++i){
			scanf("%lld",&B[n-i-1]);
			ans+=B[n-i-1]*B[n-i-1];
		}
		for(int i=0; i<n; ++i){
			A[i+n]=A[i];
			B[i+n]=0;
		}
		Convolution(A,B,2*n);
		long long mx=0;
		for(int i=n; i<2*n; ++i){
			mx=max(mx,A[i]);
		}
		printf("%lld
",ans-2*mx);
	}
	return 0;
}