jzoj4330. 几何题

Description

Input

Output
输出题目中要求的数，精确到10^-10。

Sample Input
10 5
45 70 41
9 1 43
1 68 8
70 76 7
1 19 33
71 70 53
42 54 71
11 13 30
16 63 25
30 24 34
56 61 29 7328
63 32 18 365
37 41 11 2332
36 19 43 7432
68 55 46 6338

Sample Output
6.6923868755
1.3236515510
2.2698171860
6.7830383184
5.8164492699

Data Constraint


Hint
输入输出常数巨大, 请注意优化.

题解

几(F)何(F)题(T)
因为值域很小，所以差只有-76~76共153种
把xyz压成一个153进制数s，设A[i]=∑[s[j]=i]，B[i]=∑[-s[j]=i]，把A、B移到0，卷积即可
注意会卡常

code

#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define abs(x) ((x)>0?(x):-(x))
using namespace std;

struct C{
	double x,y;
	C () {}
	C (double _x,double _y) {x=_x,y=_y;}
} a[4194304],b[4194304],A[4194304],w[4194304],W[4194304];
C operator +(C a,C b) {return C(a.x+b.x,a.y+b.y);}
C operator -(C a,C b) {return C(a.x-b.x,a.y-b.y);}
C operator *(C a,C b) {return C(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
long long ans[4194304];
int X[4194304];
int Y[4194304];
int Z[4194304];
int n,Q,i,j,k,l,N,len,x,y,z,s1,s2,s3,s4,mx,s;
double Ans;

void dft(C *a,int type)
{
	int i,j,k,l,s1=1,s2=1,S=N;
	
	fo(i,0,N-1)
	{
		j=i;
		k=0;
		fo(l,1,len)
		{
			k=k*2+(j&1);
			j>>=1;
		}
		
		A[i]=a[k];
	}
	fo(i,0,N-1)
	a[i]=A[i];
	
	l=0;
	fo(i,1,len)
	{
		s1<<=1;
		S>>=1;
		
		fo(j,0,S-1)
		{
			fo(k,0,s2-1)
			{
				C u=a[j*s1+k],v;
				if (type==1)
				v=a[j*s1+s2+k]*w[l+k];
				else
				v=a[j*s1+s2+k]*W[l+k];
				
				a[j*s1+k]=u+v;
				a[j*s1+s2+k]=u-v;
			}
		}
		
		l+=s2;
		s2<<=1;
	}
}

int main()
{
	freopen("geometry.in","r",stdin);
	freopen("geometry.out","w",stdout);
	
	l=-1;
	j=1;
	s=2;
	
	fo(i,1,22)
	{
		fo(k,0,j-1)
		{
			++l;
			w[l]=C(cos(2*M_PI*k/s),sin(2*M_PI*k/s));
			W[l].x=w[l].x;
			W[l].y=-w[l].y;
		}
		
		j<<=1;
		s<<=1;
	}
	
	scanf("%d%d",&n,&Q);
	fo(i,1,n)
	{
		scanf("%d%d%d",&x,&y,&z);
		
		++a[(x-1)*23409+(y-1)*153+(z-1)].x;
		x=77-x;
		y=77-y;
		z=77-z;
		++b[x*23409+y*153+z].x;
	}
	
	len=ceil(log(3581577)/log(2));
	N=pow(2,len);
	
	dft(a,1);
	dft(b,1);
	
	fo(i,0,N-1)
	a[i]=a[i]*b[i];
	
	dft(a,-1);
	
	fo(i,0,N-1)
	{
		ans[i]=floor(a[i].x/N+0.0001);
		
		if (ans[i])
		{
			mx=i;
			
			Z[i]=i%153-76;
			Y[i]=(i/153)%153-76;
			X[i]=(i/23409)%153-76;
		}
	}
	
	for (;Q;--Q)
	{
		scanf("%d%d%d%d",&s1,&s2,&s3,&s4);
		Ans=0;
		
		fo(i,0,mx)
		if (ans[i])
		{
			if (!X[i] && !Y[i] && !Z[i])
			continue;
			
			Ans+=abs(s1*X[i]+s2*Y[i]+s3*Z[i]+s4)/sqrt(X[i]*X[i]*X[i]*X[i]+Y[i]*Y[i]*Y[i]*Y[i]+Z[i]*Z[i]*Z[i]*Z[i])*ans[i];
		}
		
		printf("%0.6lf\n",Ans/n/(n-1));
	}
}