[JZOJ4330] 【清华集训模拟】几何题

题目

题目大意

也懒得解释题目大意了……

---

正解

正解居然是\(FFT\)？
不要看题目的那个式子这么长，也不要在那个式子上下手。
其实我们会发现，不同的\((x_i-x_j,y_i-y_j,z_i-z_j)\)并不多。
如果我们求出每个三元组的出现次数，后面的就好做了。
那怎么求呢？
祭出我们的大杀器——\(FFT\)。
考虑只有一个维怎么做。设两个多项式分别为\(A\)和\(B\)。
对于\(x_i\)，就在\(A\)的\(x_i\)这一位上的系数加一；
对于\(x_j\)，就在\(B\)的\(77-x_j\)这一位上的系数加一。
将\(A\)和\(B\)乘起来，那么\(77+x_i-x_j\)就是差\(x_i-x_j\)对应的个数。
对于三维，就将这三个数压成一维的就好了。

实际上也可以用NTT。仔细分析一下，就可以发现每个三元组的出现次数肯定是不超过\(998244353\)的。

---

正解

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cassert>
#include <cmath>
#define N 1000000
#define MX 3652264
#define mo 998244353
inline int input(){
	char ch=getchar();
	while (ch<'0' || '9'<ch)
		ch=getchar();
	int x=0;
	do{
		x=x*10+ch-'0';
		ch=getchar();
	}
	while ('0'<=ch && ch<='9');
	return x;
}
inline int my_pow(int x,int y){
	int res=1;
	for (;y;y>>=1,x=(long long)x*x%mo)
		if (y&1)
			res=(long long)res*x%mo;
	return res;
}
inline int pow4(int x){x*=x;return x*x;}
#define M (1<<22)
#define bit 22
int n;
struct DOT{
	int x,y,z;
	inline DOT rev(){return {77-x,77-y,77-z};}
} d[N];
inline int pia(const DOT &a){return (a.x*154+a.y)*154+a.z;}
int a[1<<22],b[1<<22],cnt[1<<22];
int rev[1<<22];
inline void ntt(int *a,int flag){
	for (int i=0;i<M;++i)
		if (i<rev[i])
			swap(a[i],a[rev[i]]);
	for (int i=1;i<M;i<<=1){
		int wn=my_pow(3,(mo+1)/(i<<1));
		if (flag==-1)
			wn=my_pow(wn,mo-2);
		for (int j=0;j<M;j+=i<<1){
			int wnk=1;
			for (int k=j;k<j+i;++k,wnk=(long long)wnk*wn%mo){
				int x=a[k],y=(long long)wnk*a[k+i]%mo;
				a[k]=(x+y>=mo?x+y-mo:x+y);
				a[k+i]=(x-y<0?x-y+mo:x-y);
			}
		}
	}
	if (flag==-1){
		int invm=my_pow(M,mo-2);
		for (int i=0;i<M;++i)
			a[i]=(long long)a[i]*invm%mo;
	}
}
inline void multi(int *a,int *b,int *c){
	for (int i=1;i<M;++i)
		rev[i]=rev[i>>1]>>1|(i&1)<<bit-1;
	ntt(a,1),ntt(b,1);
	for (int i=0;i<M;++i)
		c[i]=(long long)a[i]*b[i]%mo;
	ntt(c,-1);
}
DOT back[M];
int main(){
	freopen("geometry.in","r",stdin);
	freopen("geometry.out","w",stdout);
	int Q;
	scanf("%d%d",&n,&Q);
	for (int i=1;i<=n;++i)
		d[i]={input(),input(),input()};
	for (int i=1;i<=n;++i){
		a[pia(d[i])]++;
		b[pia(d[i].rev())]++;
	}
	multi(a,b,cnt);
	for (int i=0;i<MX;++i){
		int j=i;
		back[i].z=j%154-77;j/=154;
		back[i].y=j%154-77;j/=154;
		back[i].x=j-77;
//		assert(pia(back[i])==i);
	}
	while (Q--){
		int a=input(),b=input(),c=input(),d=input();
		double ans=0;
		for (int i=0;i<MX;++i)
			if (cnt[i] && (back[i].x|back[i].y|back[i].z))
				ans+=(long long)cnt[i]*abs(a*back[i].x+b*back[i].y+c*back[i].z+d)/sqrt(pow4(back[i].x)+pow4(back[i].y)+pow4(back[i].z));
		ans/=(long long)n*(n-1);
		printf("%.10lf\n",ans);
	}
	return 0;
}

---

总结

\(FFT\)和\(NTT\)真是个bug般的存在……