【hihocoder#1388】Periodic Signal NTT

题目链接:http://hihocoder.com/problemset/problem/1388?sid=974337

题目大意:找出一个$k$,使得$\sum_{i=0}^{n-1}(A_{i}-B_{i+k \quad mod \quad n})^{2}$最小

把那个式子拆开..得到:

$\sum _{i=0}^{n-1} A_{i}^{2}+\sum_{i=0}^{n-1}B_{i}^{2}-\sum_{i=0}^{n-1}2*A_{i}*B_{i+k}$

前面的两个和式是定值,所以最小化后面的值,后面这个是个卷积的形式,所以把$B$翻转硬上就好了...

然后问题就来了..FFT了之后貌似精度爆炸了..然后就写NTT..然后又炸了..发现模数规模不够..于是查表换了个巨大的模数,然后又炸了..

然后去看题解,意识到一个严重的问题..相乘爆longlong...然后又学习了TA爷的$O(1)$快速乘..

附上一张神表:折跃

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
using namespace std;
#define LL long long
inline LL read()
{
LL x=0,f=1; char ch=getchar();
while (ch<'0' || ch>'9') {if (ch=='-') f=-1; ch=getchar();}
while (ch>='0' && ch<='9') {x=x*10+ch-'0'; ch=getchar();}
return x*f;
}
#define P (LL)((29LL<<57)+1)
#define G 3
#define MAXN 800100
int T,N,len;
LL Wn[30],A[MAXN],B[MAXN],C[MAXN],wn[MAXN],a[MAXN],b[MAXN];
inline LL Mul(LL x,LL y) {return (x*y-(LL)(x/(long double)P*y+1e-3)*P+P)%P;}
inline LL Pow(LL x,LL y)
{
LL re=1;
for (LL i=y; i; i>>=1,x=Mul(x,x))
if (i&1) re=Mul(re,x);
return re;
}
inline LL Inv(LL x) {return Pow(x,P-2);}
inline void Prework()
{
len=1;
while (len<(N<<1)) len<<=1;
for (int i=0; i<N; i++) A[i]=a[i];
for (int i=0; i<N; i++) B[i]=b[N-1-i];
for (int i=N; i<len; i++) A[i]=0;
for (int i=N; i<len; i++) B[i]=0;
// for (int i=0; i<len; i++) printf("%I64d ",A[i]); puts("");
// for (int i=0; i<len; i++) printf("%I64d ",B[i]); puts("");
for (int i=0; i<=28; i++) wn[i]=Pow(G,(P-1)/(1<<i));
}
inline void Rader(LL *x)
{
for (int i=1,j=len>>1,k; i<len-1; i++)
{
if (i<j) swap(x[i],x[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
inline void DFT(LL *x,int opt)
{
Rader(x);
for (int h=2,id=0; h<=len; h<<=1)
{
LL Wn=wn[++id];
for (int i=0; i<len; i+=h)
{
LL W=1;
for (int j=i; j<i+h/2; j++)
{
LL u=x[j]%P,t=Mul(W,x[j+h/2]);
x[j]=(u+t)%P; x[j+h/2]=(u-t+P)%P;
W=Mul(W,Wn);
}
}
}
if (opt==-1)
{
for (int i=1; i<len/2; i++) swap(x[i],x[len-i]);
for (int i=0; i<len; i++)
x[i]=Mul(x[i],Inv(len));
}
}
inline void NTT(LL *A,LL *B)
{
DFT(A,1); DFT(B,1);
for (int i=0; i<len; i++) C[i]=Mul(A[i],B[i]);
DFT(C,-1);
}
int main()
{
T=read();
while (T--)
{
N=read();
LL ans=0;
for (int i=0; i<=N-1; i++) a[i]=read();
for (int i=0; i<=N-1; i++) b[i]=read();
for (int i=0; i<=N-1; i++) ans+=(LL)a[i]*a[i];
for (int i=0; i<=N-1; i++) ans+=(LL)b[i]*b[i];
Prework();
NTT(A,B);
LL mx=C[N-1];
for (int i=0; i<N-1; i++) mx=max(mx,C[i]+C[i+N]);
// for (int i=0; i<len; i++) printf("%I64d ",C[i]); puts("");
printf("%lld\n",ans-2*mx);
}
return 0;
}

  

上一篇:从程序员到CTO的Java技术路线图


下一篇:MYSQL基础03(日期函数)