JZOJ 4330 几何题

几何题

Description

给出$n$个三维空间内的点，第$i$个点的坐标是($x_i$,$y_i$,$z_i$)，现在有$q$组询问，每次询问给出$4$个数$a$,$b$,$c$,$d$，对于每次询问，求出

$$\sum_{i≠j}{{|a(x_i-x_j)+b(y_i-y_j)+c(z_i-z_j)+d|}\over{n(n-1)\sqrt{(x_i-x_j)^4+(y_i-y_j)^4+(z_i-z_j)^4}}}$$
答案保留 $6$位小数

Data Constraint

$2$≤ $n$≤$777777$
$1$≤$q$, $a$,$b$, $c$,$x_i$, $y_i$, $z_i$≤ $77$
$1$≤$d$≤ $7777$
愿数字 $7$给你带来好运。

Solution

大佬出的题数据范围就是清奇。
假设我们已经知道了$x_i-x_j$= $v$,$y_i-y_j$= $u$,$z_i-z_j$= $w$的对数，记作$A_{v,u,w}$，那么每一次询问所有可能的 $v$,$u$, $w$都枚举一次即可求出答案。

现在问题就是如何求解$A$。
我们构造这样两个多项式
$F$($x$, $y$,$z$)= $\sum_{i=1}^n x^{x_i}y^{y_i}z^{z_i}$
$G$($x$, $y$,$z$)= $\sum_{i=1}^n x^{77-x_i}y^{77-y_i}z^{77-z_i}$
显然有 $A_{v,u,w}$=[ $x^{v-77}y^{u-77}z^{w-77}$] $F(x,y,z)G(x,y,z)$
接下来只需要完成 $F(x,y,z)$和 $G(x,y,z)$的卷积即可。

考虑将三维多项式压成一维多项式。
令$F$($x$)=$\sum_{i=1}^n x^{x_i}y^{y_i*78}z^{z_i*78^2}$
$G$同理。
其实就是每一个数表示一个有三位的$78$进制数，这样便巧妙地完成压维。
实现比较卡，要卡卡常数。

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#pragma GCC optimize(3) 

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)
#define sx(o) (o)*(o)*(o)*(o)
#define ab(o) (o<0?-o:o)

using namespace std;
typedef long long ll;
typedef double db;
const ll M=84e5,N=8e5,ZD=8388608,zx=155;
const db pi=acos(-1);

struct Z{
    db x,y;
    Z(db _x=0,db _y=0){x=_x;y=_y;}
}t1[M],t2[M];

Z operator+(Z a,Z b){return Z(a.x+b.x,a.y+b.y);}
Z operator-(Z a,Z b){return Z(a.x-b.x,a.y-b.y);}
Z operator*(Z a,Z b){return Z(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

int x[M],y[M],z[M],xx,yy,zz,xy,a,b,c;
int n,q,ss,mm,x1,x2,x3,x4;
int bits[M];
ll bl[M];
db ds[zx][zx][zx];

inline void FFT_prepare()
{
//  fo(i,0,2*mm)tri[i].x=cos(i*pi/mm),tri[i].y=sin(i*pi/mm);
    fo(i,0,mm)bits[i]=(bits[i>>1]>>1)|((i&1)<<ss-1);
}

inline int read()
{
    int o=0; char ch=' ';
    for(;ch<'0'||ch>'9';ch=getchar());
    for(;ch>='0'&&ch<='9';ch=getchar())o=o*10+ch-48;
    return o;
}

inline int max(int a,int b)
{return a>b?a:b;}

void DFT(Z *a,int sig)
{
    Z tr;
    fo(i,0,mm-1)
    if(i<bits[i]){tr=a[i]; a[i]=a[bits[i]]; a[bits[i]]=tr;}
    Z u,v,yy;
    for(register int m=2,half=1;m<=mm;m<<=1,half<<=1)
        for(register int j=0,k=half;j<mm;j+=m,k+=m){
            v.x=1; v.y=0; 
            yy.x=cos(sig*pi/half);
            yy.y=sin(sig*pi/half);
            for(register int i=j,p=j+half;i<k;++i,++p,v=v*yy){
                u=v*a[p];
                a[p]=a[i]-u;
                a[i]=a[i]+u;
            }
        }
    if(sig==-1)fo(i,0,mm-1)a[i].x/=mm;
}

inline double js(int a,int b,int c)
{return sqrt(sx(a)+sx(b)+sx(c));}

inline int aabb(int o)
{return o<0?-o:o;}

int main()
{
    cin>>n>>q;
    fo(i,1,n)x[i]=read(),y[i]=read(),z[i]=read();
    fo(i,1,n)xx=max(xx,x[i]),yy=max(yy,y[i]),zz=max(zz,z[i]);
    a=(xx+1)<<1; b=(yy+1)<<1; c=(zz+1)<<1;
    --a; --b;
    xy=a*b;
    fo(i,1,n)t1[z[i]*xy+y[i]*a+x[i]].x=t1[z[i]*xy+y[i]*a+x[i]].x+1;
    fo(i,1,n)t1[(zz-z[i])*xy+(yy-y[i])*a+(xx-x[i])].y=t1[(zz-z[i])*xy+(yy-y[i])*a+(xx-x[i])].y+1;
    int zc=0;
    fo(i,1,n)zc=max(zc,z[i]*xy+y[i]*a+x[i]);
    fo(i,1,n)zc=max(zc,(zz-z[i])*xy+(yy-y[i])*a+(xx-x[i]));
    zc<<=1;
    ss=0,mm=1;
    while(mm<=zc)mm<<=1,++ss;
    FFT_prepare();
    DFT(t1,1);
    register db g1,g2,g3,g4;
    fo(i,0,mm-1){
        int j=(mm-1)&(mm-i);
        g1=(t1[i].x+t1[j].x)/2; g2=(t1[i].y-t1[j].y)/2;
        g3=(t1[i].y+t1[j].y)/2; g4=(t1[j].x-t1[i].x)/2;
        t2[i].x=g1*g3-g2*g4;
        t2[i].y=g1*g4+g2*g3;
    }
    DFT(t2,-1);
    fo(i,0,mm)bl[i]=(ll)(t2[i].x+0.3);
    fo(i3,0,zz)
    fo(i2,0,yy)
    fo(i1,0,xx)
    if(i1|i2|i3)
    ds[i3][i2][i1]=1.0/js(i3,i2,i1);
    fo(l,1,q){
        scanf("%d%d%d%d",&x1,&x2,&x3,&x4);
        register db ans=0; 
        register int i=0;
        for(register int y3=-zz;y3<=zz&&i<=mm;++y3)
        for(register int y2=-yy;y2<=yy;++y2)
        for(register int y1=-xx;y1<=xx;++y1,++i)
        if(bl[i])ans=ans+bl[i]*ds[ab(y3)][ab(y2)][ab(y1)]*aabb(x1*y1+x2*y2+x3*y3+x4);
        ans=(ans/(n-1))/n;
        printf("%.10lf\n",ans);
    }
}