Periodic Signal HihoCoder - 1388 （NTT）FFT

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

扫描二维码关注公众号，回复： 3234759 查看本文章

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

感觉这些个NTT 和FFT ，就像是树状数组一样吧，都是加快了我们的计算，至于具体的原理，不清楚

这道题，让我们求 $min$ , 通过化简，我们可以发现，其实我们只需要求出 $\sum_{i=0}^{n-1} ({A_{i}}*B_{(k+i)mod(n)})$ 的最大值即可，

然后就是板子了

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


//NTT
//begin
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);
}
//end


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]=0;
            B[i+n]=B[i];
        }

        Convolution(A,B,2*n);
        long long mx=0;
        for(int i=n-1; i<=2*n-2; ++i) {
            mx=max(mx,A[i]);
        }
        printf("%lld\n",ans-2*mx);
    }
    return 0;
}

#include<bits/stdc++.h>
using namespace std;

#define rep(i,a,b) for(int i=a;i<b;i++)
typedef long long LL;
const double pi=acos(-1);
const int maxn=(1<<20)+1;
const double eps=1e-6;

struct Com {
	double r,i;
	Com(double _r=0,double _i=0):r(_r),i(_i) {}
	Com operator+(Com& b) {
		return Com(r+b.r,i+b.i);
	}
	Com operator-(Com& b) {
		return Com(r-b.r,i-b.i);
	}
	Com operator*(Com& b) {
		return Com(r*b.r-i*b.i,i*b.r+r*b.i);
	}
} a[maxn<<2],b[maxn<<2];
int R[maxn<<2];



//FFT实现
inline void FFT(Com *x,int n,int f) {
	//Rader(x,n);
	for (int i=0; i<n; ++i) if (i<R[i]) swap(x[i],x[R[i]]);
	for (int i=1; i<n; i<<=1) {
		Com wn=Com(cos(pi/i),f*sin(pi/i));
		for (int j=0; j<n; j+=i<<1) {
			Com w=Com(1,0),t1,t2;
			for (int k=0; k<i; ++k,w=w*wn)
				t1=x[j+k],t2=x[j+k+i]*w,x[j+k]=t1+t2,x[j+k+i]=t1-t2;
		}
	}
	if(f==-1) for (int i=0; i<n; ++i) x[i].r/=n;
}

//求卷积
void Conv(Com a[],Com 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);
}

LL aa[maxn*2],bb[maxn*2];
int main() {
	int T;
	scanf("%d",&T);
	while(T--) {
		int n;
		scanf("%d",&n);
		//init 
		
		LL ans=0;
		rep(i,0,n){
			scanf("%lld",&aa[i]);
			ans+=aa[i]*aa[i];
		}
		rep(i,0,n){
			scanf("%lld",&bb[n-1-i]);
			ans+=bb[n-1-i]*bb[n-1-i]; 
		}
		rep(i,n,2*n){
			aa[i]=0,bb[i]=bb[i-n];	
		}
		//注意赋值，实数和虚数都得初始化 
		rep(i,0,2*n){
			a[i]=Com(aa[i],0);
			b[i]=Com(bb[i],0);
		}
		n*=2;//每个项式的 位数 
		
		//补0
		int len=n+n;
		int t=0,tot=1;
		for (tot=1; tot<=len; ++t,tot<<=1);
		for (int i=0; i<tot; ++i) R[i]=(R[i>>1]>>1)|(i&1)<<t-1;
		rep(i,n,tot)a[i]=Com(0,0),b[i]=Com(0,0);
			
		//卷积
		Conv(a,b,tot);

		n/=2;
	    double res=-1e18; int pos=n;
	    rep(i,n,2*n){
			if(a[i].r-res>=eps){
				res=a[i].r;
				pos=i;
			} 
		}
	
		
		LL sum=0;//int j=pos;
		rep(i,0,n){
			sum=sum+aa[i]*bb[pos-i];
		}
		ans=ans-sum*2;
		printf("%lld\n",ans);
	}
	return 0;
}

Periodic Signal HihoCoder - 1388 （NTT）FFT

猜你喜欢