几何题

题目大意

给你n个三维空间的点（坐标为整数且非负）。
q次询问，每次询问一个很神秘的式子。

做法

显然是上三维FFT。
然后需要加入一系列优化卡常。

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
//#define abs(x) (x<0?-x:x)
using namespace std;
typedef long long ll;
typedef double db;
const db pi=acos(-1);
const int maxlen=8000000+10,mo=998244353,GG=3;
int rev[maxlen],x[maxlen],y[maxlen],z[maxlen];
int sta[maxlen][3];
db si[maxlen];
int w[maxlen],a[maxlen],b[maxlen],tt[maxlen];
int i,j,k,l,r,t,n,m,len,mx,top,xx,yy,zz,ni;
db ce,ans;
int read(){
    int 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;
}
int qsm(int x,int y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%mo;
    if (y%2) t=(ll)t*x%mo;
    return t;
}
void prepare(){
    fo(i,0,len-1){
        int p=0;
        for(int j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
        rev[i]=p;
    }
    w[0]=1;
    w[1]=qsm(GG,(mo-1)/len);
    fo(i,2,len) w[i]=(ll)w[i-1]*w[1]%mo;
    ni=qsm(len,mo-2);
}
void DFT(int sig){
    int i;
    fo(i,0,len-1) tt[rev[i]]=a[i];
    for(register int m=2;m<=len;m*=2){
        /*int half=m/2,bei=len/m;
        fo(i,0,half-1){
            int wi=sig>0?w[i*bei]:w[len-i*bei];
            for(int j=i;j<len;j+=m){
                int u=tt[j],v=(ll)tt[j+half]*wi%mo;
                tt[j]=(u+v)%mo;
                tt[j+half]=(u-v)%mo;
            }
        }*/
        register int half=m/2;
        register int wi=sig>0?w[len/m]:w[len-len/m];
        for(i=0;i<len;i+=m){
            int o=1;
            for(register int j=i;j<i+half;j++,o=(ll)o*wi%mo){
                int v=(ll)o*tt[j+half]%mo;
                tt[j+half]=(tt[j]-v+mo)%mo;
                //tt[j]=(tt[j]+v)%mo;
                tt[j]=tt[j]+v>=mo?tt[j]+v-mo:tt[j]+v;
            }
        }
    }
    if (sig==-1)
        fo(i,0,len-1) tt[i]=(ll)tt[i]*ni%mo;
    fo(i,0,len-1) a[i]=tt[i];
}
void NTT(){
    int i;
    DFT(1);
    fo(i,0,len-1) swap(a[i],b[i]);
    DFT(1);
    fo(i,0,len-1) a[i]=(ll)a[i]*b[i]%mo;
    DFT(-1);
}
void work(int a,int b,int c,int d){
    swap(a,d);
    swap(b,c);
    int i,t,x,y,z;
    ans=0;
    fo(i,1,top){
        /*x=sta[i][0];
        y=sta[i][1];
        z=sta[i][2];*/
        //t=sta[i][3];
        //ans+=(db)abs(a*x+b*y+c*z+d)/si[i]/**t*/;
        ans+=(db)abs(a*sta[i][0]+b*sta[i][1]+c*sta[i][2]+d)*si[i]/**t*/;
    }
    ans=ans/n/(n-1);
    printf("%.15lf\n",ans);
}
int main(){
    freopen("geometry15.in","r",stdin);freopen("geometry.out","w",stdout);
    n=read();m=read();
    fo(i,1,n){
        x[i]=read();y[i]=read();z[i]=read();
        mx=max(mx,x[i]);
        mx=max(mx,y[i]);
        mx=max(mx,z[i]);
    }
    //mx=77;
    len=1;
    r=2*mx+1;
    l=2*mx*r*r+2*mx*r+2*mx;
    while (len<=l) len*=2;
    ce=log(len)/log(2);
    prepare();
    fo(i,1,n){
        t=x[i]*r*r+y[i]*r+z[i];
        a[t]++;
        t=(mx-x[i])*r*r+(mx-y[i])*r+(mx-z[i]);
        b[t]++;
    }
    NTT();
    //return 0;
    fo(i,0,len-1){
        l=a[i];
        (l+=mo)%=mo;
        if (!l) continue;
        t=i;
        xx=t/(r*r)-mx;
        t%=(r*r);
        yy=t/r-mx;
        t%=r;
        zz=t-mx;
        if (xx==0&&yy==0&&zz==0) continue;
        top++;
        sta[top][0]=xx;
        sta[top][1]=yy;
        sta[top][2]=zz;
        //sta[top][3]=l;
        si[top]=sqrt(xx*xx*xx*xx+yy*yy*yy*yy+zz*zz*zz*zz)/l;
        si[top]=1/si[top];
    }
    //return 0;
    while (m--) work(read(),read(),read(),read());
}