Anarchy的解题报告

题目大意：

假设高斯定理在m维空间成立

已知m维空间所有整点电荷 aj

给出m维空间下x,y两点距离公式



以及x点在y点引发的电势公式





降序输出求前100个点的电势

T≤51≤m≤18n≤3∗105

时限7s

考试的时候看这道题的时候 感觉这道题根本不可做啊？？？

部分分很友好啊？除了一个5分的点 其他我都不会。。

看完题解之后整个人都不好了。。

然后我用了一个下午来调试 压常

其实 上面那个式子就是卷积的形式 只不过 这个是m维卷积。。。（做毛线啊

好吧 m维卷积大家都会喜闻乐见的 FWT ？

平常的FWT一般都是每一维大小为2 的情况 当大于2的时候我们可以直接上DFT

然后我们会发现这个循环卷积也不好处理

怎么办呢？

将DFT的式子拆成卷积形式 然后FFT优化卷积之后即可得出

即Bluestein′sAlgorithm

时间复杂度O(nlog2n)

注意常数很大 所以不得不不预处理单位根 以及一些黑科技

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
#define ld double
ld Sin[20],Cos[20];
struct com
{
ld a,b;
com(){}
com(ld a,ld b):a(a),b(b){}
friend com operator +(com a,com b){return com(a.a+b.a,a.b+b.b);}
friend com operator -(com a,com b){return com(a.a-b.a,a.b-b.b);}
friend com operator *(com a,com b){return com(a.a*b.a-a.b*b.b,a.b*b.a+a.a*b.b);}
friend com operator /(com a,ld b){return com(a.a/b,a.b/b);}
friend com operator *(com a,ld b){return com(a.a*b,a.b*b);}

};
const
ld pi=acos(-1);
int rev[2200005];
void FFT(com *a,int n,int fl)
{
int Len=1<<n;
for(int i=0;i<Len;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<n-1);
for(int i=0;i<Len;i++)if(rev[i]>i)swap(a[rev[i]],a[i]);
for(int p=0,i=1;i<Len;i*=2,p++)
{
com w=com(Cos[p],fl*Sin[p]);
for(int j=0;j<Len;j+=2*i)
{
com w0=com(1,0);
for(int k=0;k<i;k++)
{
com x=a[j+k],y=w0*a[i+j+k];
a[j+k]=x+y;
a[j+k+i]=x-y;
w0=w0*w;
}
}
}
if(fl==-1)
for(int i=0;i<Len;i++)a[i]=a[i]/Len;
}

com c[2200005],b[2200005];
int C;
ld SIN[2200005],COS[2200005];
void DFT(com* a,int d,int n,int fl)
{
int N=1;
for(;(1<<N)<(3*n);N++);
int P=C/2/n,Nf=2*n;
for(int i=0;i<2*n;i++)
b[i]=com(COS[i*1ll*i%Nf*P],fl*SIN[i*1ll*i%Nf*P]);
for(int i=0;i<n;i++)
c[n-i]=com(COS[i*1ll*i%Nf*P],-fl*SIN[i*1ll*i%Nf*P])*a[i*d];
if(n<=50)
{
for(int i=n;i<2*n;i++)
{
int pos=i-n;
pos*=d;
a[pos]=com(0,0);
for(int j=0;j<=i;j++)
a[pos]=a[pos]+b[i-j]*c[j];
}
}
else
{
for(int i=0;i<N;i++)
Cos[i]=cos(pi/(1<<i)),Sin[i]=sin(pi/(1<<i));
FFT(b,N,1);
FFT(c,N,1);
for(int i=0;i<(1<<N);i++)b[i]=b[i]*c[i];
FFT(b,N,-1);
for(int i=0;i<n;i++)a[i*d]=b[i+n];
}
for(int i=0;i<n;i++)
a[i*d]=a[i*d]*com(COS[i*1ll*i%Nf*P],-fl*SIN[i*1ll*i%Nf*P]);
if(fl==-1)
for(int i=0;i<n;i++)a[i*d]=a[i*d]/n;
for(int i=0;i<(1<<N);i++)b[i]=c[i]=com(0,0);
}
int s[2200005];
void FWT(com* a,int n,int m,int fl)
{
int len=1;
for(int i=0;i<m;len*=s[i++])
for(int j=0;j<n;j+=len*s[i])
for(int k=0;k<len;k++)
DFT(a+j+k,len,s[i],fl);
}
com Da[2200005];
int no[2200005];
bool cmp(int a,int b)
{
return Da[a].a>Da[b].a;
}
int main()
{
freopen("anarchy.in","r",stdin);
freopen("anarchy.out","w",stdout);
int T;
scanf("%d",&T);
while(T--)
{
int m,n=1;
scanf("%d",&m);
for(int i=0;i<m;i++)scanf("%d",s+i),n*=s[i];C=n;
for(int i=0;i<n;i++)
{
int t;
scanf("%d",&t);
Da[i].a=t;
Da[i].b=1;
for(int t=i,j=0;j<m;t/=s[j++])
Da[i].b*=t%s[j]+1;
Da[i].b=powl(Da[i].b,2.0/m);
Da[i].b/=2.*m;
Da[i]=com(Da[i].a+Da[i].b,Da[i].a-Da[i].b);
}
C*=2;
for(int i=0;i<2*n;i++)
SIN[i]=sin(2*pi/2/n*i),COS[i]=cos(2*i*pi/2/n);
FWT(Da,n,m,1);
for(int i=0;i<n;i++)Da[i]=Da[i]*Da[i];
FWT(Da,n,m,-1);

for(int i=0;i<n;i++)no[i]=i;
sort(no,no+n,cmp);
for(int i=0;i<min(n,100);i++)
printf("%.9f ",Da[no[i]].a/4);
puts("");

}
return 0;
}