TopCoder SRM 606 Div1 Hard Subcube

問題

空間上の点の集合がsubcubeであるとは、それらの点が、ある立方体の頂点の部分集合になっていることを言う。


n個の点が与えられる。
この点の部分集合のうち、subcubeであるようなものの個数を求めよ。

制約条件

n≦50
座標の絶対値≦100万

方針

全生成して間に合うらしい・・
4点が乗る立方体は高々2通りしかない。5点が乗る立方体は1通りしかない。
したがってそんなに組み合わせの個数が多くならない。


自分はこれに途中まで気づかず、3点に乗る立方体が少ないので、
それを全生成して、立方体からそれに乗る点を数えて、コンビネーションから場合の数を計算する、とかやってやばいことになっていた。


1点、2点からなる集合は全てsubcube
3点は、三角形の辺の長さの比が1:1:1または1:√2:√3または1:1:√2になる場合だけがsubcube
で、そこから立方体は2通りに定まるので、それを全生成。


4点以上のsubcubeは立方体のほうから数える。
立方体ごとに乗っている点の個数を数える。
C[cnt][4] + C[cnt][5] + C[cnt][6] + C[cnt][7] + C[cnt][8]が組み合わせの数。


でも実は4点は二度数えられている可能性があって、
その重複は個数が少なかったのでsetに入れたりしてがんばって取り除いた。


でもそれだったら最初から全部生成して間に合うよねっていう悲しいことになった。
立方体の生成とかやる必要なかった……

typedef vector<double> P;
P operator+(const P&a, const P&b){
	P c(3);
	rep(i, 3) c[i] = a[i] + b[i];
	return c;
}
P operator-(const P&a, const P&b){
	P c(3);
	rep(i, 3) c[i] = a[i] - b[i];
	return c;
}
P operator*(double a, const P&b){
	P c(3);
	rep(i, 3) c[i] = a * b[i];
	return c;
}
P op(const P&a, const P&b){
	P c(3);
	c[0] = a[2] * b[1] - a[1] * b[2];
	c[1] = a[0] * b[2] - a[2] * b[0];
	c[2] = a[1] * b[0] - a[0] * b[1];
	return c;
}
double norm(const P&a){
	double r = 0;
	rep(i, 3) r += a[i] * a[i];
	return r;
}
namespace std{
bool operator<(const P &a, const P &b){
	rep(i, 3){
		if(a[i] + EPS < b[i]) return 1;
		if(a[i] > b[i] + EPS) return 0;
	}
	return 0;
}}

struct C{
	P p[8];
	bool operator<(const C &o)const{
		rep(i, 8) rep(j, 3){
			if(p[i][j] + EPS < o.p[i][j]) return 1;
			if(p[i][j] > o.p[i][j] + EPS) return 0;
		}
		return 0;
	}
};


int n;
vi x, y, z;
set<C> cubes;
ll com[60][60];

double dist(int a, int b){
	double r = 0;
	r += 1.0 * (x[a] - x[b]) * (x[a] - x[b]);
	r += 1.0 * (y[a] - y[b]) * (y[a] - y[b]);
	r += 1.0 * (z[a] - z[b]) * (z[a] - z[b]);
	return r;
}
int subcube3(int a, int b, int c){
	double d[] = {dist(a, b), dist(b, c), dist(c, a)};
	sort(d, d + 3);
	if(abs(d[1] - d[0]) < EPS && abs(d[2] - d[0]) < EPS) return 1;
	if(abs(d[1] - d[0]) < EPS && abs(d[2] - d[0] * 2) < EPS) return 2;
	if(abs(d[1] - d[0] * 2) < EPS && abs(d[2] - d[0] * 3) < EPS) return 3;
	return 0;
}
void ins1(const P&a, const P& b, const P& c){
 	P g = a + b + c;
 	g = (1.0 / 3) * g;
	
	P d = op(b - a, c - a);
	d = (sqrt(norm(a - b)) / sqrt(6.0) / sqrt(norm(d))) * d;
	
	vector<P> ps;
	ps.pb(a); ps.pb(b); ps.pb(c);
	g = g + d;
	
	ps.pb(g); ps.pb(a + b - g); ps.pb(a + c - g); ps.pb(b + c - g);
	ps.pb(a + (c - g) + (b - g));
	
	sort(all(ps));
	C cu;
	rep(i, 8) cu.p[i] = ps[i];
	cubes.insert(cu);
}
void ins2(P a, P b, P c){
	if(abs(norm(b - a) - norm(c - a)) < EPS) swap(a, c);
	else if(abs(norm(a - b) - norm(c - b)) < EPS) swap(b, c);
	
	P d = op(a - c, b - c);
	d = (sqrt(norm(a - c)) / sqrt(norm(d))) * d;
	
	P m = 0.5 * (a + b);
	
	vector<P> ps;
	ps.pb(a); ps.pb(b); ps.pb(c);
	ps.pb(a + d); ps.pb(b + d); ps.pb(c + d);
	ps.pb(c + 2.0 * (m - c)); ps.pb(ps.back() + d);
	
	sort(all(ps));
	C cu;
	rep(i, 8) cu.p[i] = ps[i];
	cubes.insert(cu);
}
void ins3(P a, P b, P c){
	vector<P> ps;
	ps.pb(a); ps.pb(b); ps.pb(c);
	int o[3] = {0, 1, 2};
	do{
		bool ok = 1;
		if(abs(norm(ps[o[0]] - ps[o[1]]) * 2 - norm(ps[o[1]] - ps[o[2]])) > EPS) ok = 0;
		if(abs(norm(ps[o[0]] - ps[o[1]]) * 3 - norm(ps[o[0]] - ps[o[2]])) > EPS) ok = 0;
		if(ok) break;
	}while(next_permutation(o, o + 3));
	
	a = ps[o[0]]; b = ps[o[1]]; c = ps[o[2]];
	ps.pb(a + c - b);
	P d = op(a - b, c - b);
	d = (sqrt(norm(a - b)) / sqrt(norm(d)) / sqrt(2.0)) * d;
	ps.pb(a + 0.5 * (c - b) + d); ps.pb(a + 0.5 * (c - b) - d);
	ps.pb(0.5 * (b + c) + d); ps.pb(0.5 * (b + c) - d);
	sort(all(ps));
	
	C cu;
	rep(i, 8) cu.p[i] = ps[i];
	cubes.insert(cu);
}

void pv(ll a){
	rep(i, 60) cerr<<(a>>i&1)<<(i==59?"\n":"");
}

class Subcube {
	public:
	long long getCount(vector <int> X, vector <int> Y, vector <int> Z) {
		x = X; y = Y; z = Z;
		n = x.size();
		cubes.clear();
		
		rep(i, 60) rep(j, i+1) com[i][j] = j==0||j==i ? 1 : com[i-1][j] + com[i-1][j-1];
		
		ll ans = n + com[n][2];
		
		rep(i, n) rep(j, i) rep(k, j){
			int type = subcube3(i, j, k);
			if(type == 0) continue;
			
			P a(3), b(3), c(3);
			a[0] = X[i]; a[1] = Y[i]; a[2] = Z[i];
			b[0] = X[j]; b[1] = Y[j]; b[2] = Z[j];
			c[0] = X[k]; c[1] = Y[k]; c[2] = Z[k];
			
			if(type == 1) ins1(a, b, c);
			if(type == 2) ins2(a, b, c);
			if(type == 3) ins3(a, b, c);
			ans++;
		}
		
		vector<ll> v;
		map<ll, set<int> > dupfour;
		
		each(i, cubes){
			ll bit = 0;
			rep(j, n){
				bool ok = 0;
				rep(k, 8)
				if(abs(i->p[k][0] - X[j]) < EPS &&
					abs(i->p[k][1] - Y[j]) < EPS &&
					abs(i->p[k][2] - Z[j]) < EPS){
					ok = 1;
					break;
				}
				if(ok) bit |= 1ll << j;
			}
			int cnt = __builtin_popcountll(bit);
			if(cnt > 3) v.pb(bit);
			for(int j = 4; j <= cnt; j++) ans += com[cnt][j];
		}
		rep(i, v.size()) rep(j, i){
			ll b = v[i] & v[j];
			/*
			assert(__builtin_popcountll(b) <= 4);
			if(__builtin_popcountll(b) > 4){
				pv(v[i]);
				pv(v[j]); cerr<<endl;
			}
			*/
			if(__builtin_popcountll(b) == 4){
				dupfour[b].insert(i);
				dupfour[b].insert(j);
			}
		}
		each(i, dupfour) ans -= i->second.size() - 1;
		return ans;
	}
};