【JZOJ 4330】【清华集训模拟】几何题

Description

Solution

这题的复杂度是$O(10^6\log(10^6))$的FFT（害怕.jpg）
预处理cnt[i][j][k]表示x,y,z的差为i,j,k的有多少对，
首先，$x_i-x_j$可以变成$x_i+(mx-x_j)$，（mx为最大$x_i$）
那么每位都是非负数了，
现在要做3维的多项式乘法，
考虑用2mx进制存储这三个数，压在一起变一个数，这样就把3维变成1维了，直接用FFT
剩下的就是卡常了，看我的代码吧，

复杂度：$O((77*2)^3*\log((77*2)^3))$

Code

#include <cstdio>
#include <algorithm>
#include <cmath>
#define fo(i,a,b) for(register int i=a;i<=b;++i)
#define min(q,w) ((q)>(w)?(w):(q))
#define max(q,w) ((q)<(w)?(w):(q))
#define pre(q) ((q)*(q))
#define JS(q,w,e) ((q)*pre(mx+1)*4+(w)*(mx+1)*2+(e))
using namespace std;
typedef double db;
typedef long long LL;
const int M=(1<<22)+5;
const db PI=acos(-1);
int read(int &n)
{
    char ch=' ';int q=0,w=1;
    for(;(ch!='-')&&((ch<'0')||(ch>'9'));ch=getchar());
    if(ch=='-')w=-1,ch=getchar();
    for(;ch>='0' && ch<='9';ch=getchar())q=q*10+ch-48;n=q*w;return n;
}
int m1,n,mx;
db ans;
struct Z
{
    db x,y;
    Z(db _x=0,db _y=0){x=_x;y=_y;}
    friend Z operator +(Z q,Z w){return Z(q.x+w.x,q.y+w.y);}
    friend Z operator -(Z q,Z w){return Z(q.x-w.x,q.y-w.y);}
    friend Z operator *(Z q,Z w){return Z(q.x*w.x-q.y*w.y,q.x*w.y+q.y*w.x);}
}c1[M],c[M],ansf[M];
int b1[M],b2[M];
db Ans[M];
void DFT(Z *a,int n,int K)
{
    fo(i,0,n-1)
    {
        int q=0;
        for(int j=i,k=n-1;k;k>>=1,j>>=1)q=(q<<1)+(j&1);
        c[q]=a[i];
    }
    for(int I=2;I<=n;I<<=1)
    {
        register int mid=I>>1;
        Z W(cos(K*PI/mid),sin(K*PI/mid));
        for(int j=0;j<n;j+=I)
        {
            Z w(1.0,0);
            fo(i,j,j+mid-1)
            {
                Z t=c[i+mid]*w;
                c[i+mid]=c[i]-t;
                c[i]=c[i]+t;
                w=w*W;
            }
        }
    }
    if(K<0)fo(i,0,n-1)c[i].x=c[i].x/n;
}
void FFT()
{
    int m=(1<<22);
    fo(i,0,m)c1[i]=Z(b1[i],b2[i]);
    DFT(c1,m,1);
    db ax,ay,bx,by;
    fo(i,0,m-1)
    {
        int j=(m-i)&(m-1);
        ax=(c[i].x+c[j].x)/2,ay=(c[i].y-c[j].y)/2;
        bx=(c[i].y+c[j].y)/2,by=(c[j].x-c[i].x)/2;
        c1[i].x=ax*bx-ay*by;
        c1[i].y=ax*by+ay*bx;
    }
    DFT(c1,m,-1);
}
int main()
{
    freopen("geometry.in","r",stdin);
    freopen("geometry.out","w",stdout);
    int q,w,e;
    read(n),read(m1);
    mx=77;
    fo(i,1,n)
    {
        read(q),read(w),read(e);
        ++b1[JS(q,w,e)];
        ++b2[JS(mx-q,mx-w,mx-e)];
    }
    FFT();
    fo(i,0,2*mx)
    {
        register int x=i-mx;
        fo(j,0,2*mx)
        {
            register int y=j-mx;
            fo(k,0,mx*2)
            {
                if(pre(x)+pre(y)+pre(k-mx)!=0)
        Ans[JS(i,j,k)]=(1.0/sqrt(pre(pre(x))+pre(pre(y))+pre(pre(k-mx))))*((LL)(c[JS(i,j,k)].x+0.4));
            }
        }
    }
    fo(I,1,m1)
    {
        register int t1,t2,t3,t4;
        t1=read(q),t2=read(q),t3=read(q),t4=read(q);
        ans=0;
        fo(i,0,2*mx)
        {
            register int x=i-mx;
            fo(j,0,2*mx)
            {
                register int y=j-mx;
                fo(k,0,mx*2)
                {
                    ans+=abs(t1*x+t2*y+t3*(k-mx)+t4)*Ans[JS(i,j,k)];
                }
            }
        }
        printf("%.10lf\n",ans/(n-1.0)/n);
    }
    return 0;
}