2016ICPC北京站网络预选赛F hihoCoder1388 Periodic Signal,FFT|NTT+CRT

// author:latstars
// time:26/092016 19:52
// Problem:hihoCoder #1388 : Periodic Signal
// Class:FFT plus special judge 
// or NTT + CRT
//  题意:求最小的 sum(A[i]-B[i+k% n])^2  for i=0 to n-1 for k=0 to n-1 
//  思路:展开 得到sum (A[i]^2+B[i+k %n]^2-2*A[i]*B[i+k %n])  for i=0 to n-1
//  前面的A[i]^2 和 B[i+k %n]^2 的sum 是一个定值
//  那么只要求sum(A[i]*B[i+k%n]) for i=0 to n-1的最大值就可以
// 是一个卷积
// 构造多项式乘法
// A[0] A[1] A[2] ...... A[n-1] 0 0 0 ...... 0 0 0
// B[n-1] B[n-2] ......  B[0]  B[n-1] B[n-2] ...... B[2] B[1] B[0] 
// 所以 呢 C[n-1] =A[0]B[0]+A[1]B[1] ....
// C[n]=A[0]B[n-1]+A[1]B[0]+...
// C[2*n-2]=A[0]B[1]+A[1]B[2]+...
// 但是会出现精度误差,那么可以这样找出最大的卷积对应的位置
// 然后再算一下原序列

// 另外一种方法是用NTT 但是不能使用p为long long 的模数,因为会暴long long 
// 所以可以选择多个不同的模数,然后使用中国剩余定理还原原来的数

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<complex>
#include<vector>
#include<algorithm>
using namespace std;
const int maxn=(1<<21)+1000;
const int maxm=(1<<18)+10;
const double PI=acos(-1.0);
typedef complex<long double> Complex;
void FFT(Complex *  a,int n,int inv)
{
    for(int i=0,j=0;i<n;i++){
        if(j>i) swap(a[i],a[j]);
        int k=n;
        while(j&(k>>=1))    j&=~k;
        j|=k;
    }
    for(int s=1;s<n;s<<=1){//s=2^0 ,s wei chuli de DFT yuanshu de yiban
        double ap=PI/s*inv;//ap=2*PI/(s*2)*inv
        Complex wn=Complex(1.0,0.0);
        Complex w=Complex(cos(ap),sin(ap));//w=
        for(int k=0;k<s;k++){               
            // cout<<wn<<endl;
            // w=Complex(cos(ap*k),sin(ap*k));
            for(int e=k;e<n;e+=s<<1){
                int o=e+s;
                Complex u=wn*a[o];
                a[o]=a[e]-u;
                a[e]+=u;
            }
            wn=wn*w;
        }
    }
    if(inv==-1){
        for(int i=0;i<n;i++) a[i]/=n;
    }
}
typedef long long LL;
LL A[maxm],B[maxm];
Complex C[maxn];
Complex D[maxn];
int T;
int n;
void solve(void)
{

}

int main()
{
    scanf("%d",&T);
    while(T--){
        LL ans=0;
        memset(A,0,sizeof(A));
        memset(B,0,sizeof(B));
        scanf("%d",&n);
        for(int i=0;i<n;i++){
            scanf("%lld",&A[i]);
            // A[i]=1<<20;
            ans=ans+A[i]*A[i];
            // cout<<ans<<endl;
        }
        for(int i=0;i<n;i++){
            scanf("%lld",&B[n-i-1]);
            // B[n-i-1]=1<<20;
            ans=ans+B[n-i-1]*B[n-i-1];
        }

        for(int i=n;i<2*n;i++){
            A[i]=0;
            B[i]=B[i-n];
        }
        n=n*2;
        int len=1;
        while(len<=n)   len<<=1;
        len<<=1;
        // cout<<len<<endl;
        // cout<<ans<<endl;
        for(int i=0;i<n;i++){
            C[i]=Complex(A[i],0.0);
            D[i]=Complex(B[i],0.0);
        }
        for(int i=n;i<len;i++){
            C[i]=D[i]=Complex(0.0,0.0);
        }
        FFT(C,len,1);
        FFT(D,len,1);
        for(int i=0;i<len;i++){
            C[i]=C[i]*D[i];
        }
        FFT(C,len,-1);
        LL tmp=0;
        int pos=0;
        for(int i=n/2;i<n;i++){//
            if(tmp<=((C[i].real()+0.1))){
                tmp=(C[i].real()+0.1);
                pos=i;
            }
            // cout<<i<<" "<<(LL)(C[i].real()+0.1)<<endl;
        }
        tmp=0;
        for(int i=0;i<n/2;i++){
               tmp+=A[i]*B[pos-i];
        }
        printf("%lld\n",ans-2*tmp);
    }

    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值