Anarchy的解题报告

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/liutian429073576/article/details/61198405

题目大意:
假设高斯定理在m维空间成立
已知m维空间所有整点电荷 aj
给出m维空间下x,y两点距离公式
这里写图片描述
以及x点在y点引发的电势公式
这里写图片描述

这里写图片描述
降序输出求前100个点的电势

T51m18n3105
时限7s

考试的时候看这道题的时候 感觉这道题根本不可做啊???
部分分很友好啊?除了一个5分的点 其他我都不会。。

看完题解之后整个人都不好了。。
然后我用了一个下午来调试 压常

其实 上面那个式子就是卷积的形式 只不过 这个是m维卷积。。。(做毛线啊
好吧 m维卷积大家都会喜闻乐见的 FWT ?
平常的FWT一般都是每一维大小为2 的情况 当大于2的时候我们可以直接上DFT
然后我们会发现这个循环卷积也不好处理
怎么办呢?

将DFT的式子拆成卷积形式 然后FFT优化卷积之后即可得出
BluesteinsAlgorithm
时间复杂度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;
}


阅读更多

没有更多推荐了,返回首页