Loj 2407 气氛/Luft

Posted on | Edited on | Views:

高维多面体的体积计算 + 高斯消元的一个 trick.

如果只有 $n$ 个 $n-1​$ 维的点,那么凸包其实就是它们围成的高维多面体.

它的体积就是选一个起点, 将 $n-1$ 个 $n-1$ 维向量形成的矩阵的行列式绝对值除掉 $(n-1)!​$ .

现在有 $n+1$ 个点,求它们形成的凸包体积.

有一个通用的结论,这 $n+1$ 个 $n-1$ 维点的凸包体积等于任选出 $n$ 个点的高维多面体体积之和 $/2​$ .

枚举哪一个点没选,用高斯消元求出其余 $n​$ 个点形成的高维多面体体积即可.

注意不能直接在模意义下计算,否则只能算出行列式在模意义下的值,但无法判定符号,得出绝对值.

所以必须要用真实值消元,为了避免 double 的精度误差,可以将真实值和模意义下的值分别维护.

最后根据真实值的符号决定是否对模意义下的值取相反数.

时间复杂度 $O(t\cdot n^4)$ .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
//%std
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out = 0, fh = 1;
	char jp = getchar();
	while ((jp > '9' || jp < '0') && jp != '-')
		jp = getchar();
	if (jp == '-')
		fh = -1, jp = getchar();
	while (jp >= '0' && jp <= '9')
		out = out * 10 + jp - '0', jp = getchar();
	return out * fh;
}
void print(int x)
{
	if (x >= 10)
		print(x / 10);
	putchar('0' + x % 10);
}
void write(int x, char c)
{
	if (x < 0)
		putchar('-'), x = -x;
	print(x);
	putchar(c);
}
const int P = 1e9 + 7;
int add(int a, int b)
{
	return a + b >= P ? a + b - P : a + b;
}
void inc(int &a, int b)
{
	a = add(a, b);
}
int mul(int a, int b)
{
	return 1LL * a * b % P;
}
int fpow(int a, int b)
{
	int res = 1;
	while (b)
	{
		if (b & 1)
			res = mul(res, a);
		a = mul(a, a);
		b >>= 1;
	}
	return res;
}
const int N = 40;
int n, p[N][N], b[N][N];
long double a[N][N];
int Guass()
{
	long double res = 1;
	int prod = 1, sgn = 1;
	for (int i = 1; i < n; ++i)
	{
		for (int j = i; j < n; ++j)
			if (b[j][i])
			{
				if (j != i)
				{
					sgn *= -1, res *= -1;
					swap(b[i], b[j]);
					swap(a[i], a[j]);
				}
				break;
			}
		if (!b[i][i])
			return 0;
		res *= a[i][i], prod = mul(prod, b[i][i]);
		int inv = fpow(b[i][i], P - 2);
		for (int j = i + 1; j < n; ++j)
			{
				long double t = -a[j][i] / a[i][i];
				int r = mul(P - b[j][i], inv);
				for (int k = 1; k < n; ++k)
					a[j][k] += t * a[i][k], inc(b[j][k], mul(r, b[i][k]));
			}
	}
	if (sgn == -1)
		prod = add(0, P - prod);
	if (res < 0)
		return add(0, P - prod);
	return prod;
}
void solve()
{
	n = read();
	for (int i = 1; i <= n + 1; ++i)
		for (int j = 1; j < n; ++j)
			p[i][j] = read();
	int ans = 0;
	for (int ban = 1; ban <= n + 1; ++ban)
	{
		int x = ban == 1 ? 2 : 1;
		int t = 0;
		for (int i = 1; i <= n + 1; ++i)
			if (i != ban && i != x)
			{
				++t;
				for (int j = 1; j < n; ++j)
                {	
                    a[t][j] = p[i][j] - p[x][j];
                    b[t][j] = add(p[i][j], P - p[x][j]);
                }
			}
		inc(ans, Guass());
	}
	ans = mul(ans, (P + 1) >> 1);
	write(ans, '\n');
}
int main()
{
	int T = read();
	while (T--)
		solve();
	return 0;
}