专栏文章

题解：P7429 [THUPC 2017] 气氛

EEdward2019
P7429题解参与者 1已保存评论 0

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
0 条
当前快照
1 份
快照标识符
@mind9ccg
此快照首次捕获于
2025/12/02 00:31
3 个月前
此快照最后确认于
2025/12/02 00:31
3 个月前
查看原文
提供一个 O(Tn3)O(Tn^3) 的爆标做法。
首先结论是 n−1n−1 维空间中 n+1n+1 个点构成的凸包的体积等于所有选 nn 个点构成的凸包体积之和的一半，其他题解对此已有详细说明。
然后考虑 n−1n−1 维空间中 nn 个点构成的凸包体积怎么算。先随便选出一个点，然后用其他 n−1n-1 个点到这个点的向量求行列式的绝对值即为体积。
首先根据 OEIS A003432，3535 阶所有元素全为 0,10,1 的行列式的上界是 2×3362\times 3^{36}，证明我还不会。这个说明我们使用 C++ 浮点数类型在本题数据范围下是不会爆精度的。
我们可以把求体积过程中求向量的“减”改为“异或”，这样就是全 0,10,1 行列式，而且也是对的。考虑这就是把空间转了一下然后把那个点转到原点，而我们只需要知道有向体积的绝对值，所以答案不变。
我们先 O(n3)O(n^3) 直接计算第 2,3,...n+12,3,...n+1 个点组成的凸包体积，然后把除了第一个向量以外的向量全部“异或”上第一个向量。
考虑下面这组数据：
CPP
1
4
1 0 1 
0 0 1 
1 1 1 
0 1 1 
0 1 0
“异或”上第一个向量后得到的四个向量是：
CPP
1 0 0
0 1 0
1 1 0
1 1 1
我们要求的就是分别去掉每一行后的行列式之和。这个形式有点像代数余子式，我们考虑给这个矩阵左边添上一列：
CPP
0 1 0 0
0 0 1 0
0 1 1 0
0 1 1 1
求一个方阵的所有代数余子式是可以做到 O(n3)O(n^3) 的（首先你需要了解伴随矩阵的概念，这里记为 adj⁡(A)\operatorname{adj}(A)），过程如下：
  • 如果 rank⁡(A)<n−1\operatorname{rank}(A)<n-1，那么显然 adj⁡(A)=0\operatorname{adj}(A)=0。
  • 如果 rank⁡(A)=n\operatorname{rank}(A)=n，那么直接使用 adj⁡(A)=det⁡(A)A−1\operatorname{adj}(A)=\det(A)A^{-1} 求出伴随矩阵即可（不过这种情况显然本题中不会出现了）。
接下来就是 rank⁡(A)=n−1\operatorname{rank}(A)=n-1 的情况了。
首先当 rank⁡(A)=n−1\operatorname{rank}(A)=n-1 时 rank⁡(adj⁡(A))=1\operatorname{rank}(\operatorname{adj}(A))=1，证明考虑，首先显然 rank⁡(adj⁡(A))≥1\operatorname{rank}(\operatorname{adj}(A))\ge1，并且有 Aadj⁡(A)=0A\operatorname{adj}(A)=0。然后有矩阵秩经典性质：对于任意矩阵 Am×nA_{m\times n} 和 Bn×sB_{n\times s}，AB=0AB=0 时 rank⁡(A)+rank⁡(B)≤n\operatorname{rank}(A)+\operatorname{rank}(B)\le n，所以 rank⁡(adj⁡(A))≤1\operatorname{rank}(\operatorname{adj}(A))\le1。
所以 adj⁡(A)\operatorname{adj}(A) 可以被表示为 uvTuv^T 的形式，其中 u,vu,v 是 nn 阶列向量。
因为 Aadj⁡(A)=adj⁡(A)A=0A\operatorname{adj}(A)=\operatorname{adj}(A)A=0，所以 u,vu,v 分别在 AA 的零空间、左零空间当中，即 Au=ATv=0Au=A^Tv=0。我们只需要从对应的解空间中分别解出任意一组非零解 u′,v′u',v'，那么容易发现 uvuv 是 u′v′u'v' 的常数倍。
在本题中已经在左边添上了全零的一列，所以可以直接取 u′=(1,0,0,...,0)u'=(1,0,0,...,0)。由于本题中我们最后只需要知道绝对值，所以可以省去很多关于 −1-1 的若干次方的讨论，那么容易发现我们最终只需要求 v′v'，然后取一个 v′v' 中的非零位置用朴素方法求一下那个位置对应的子式再除以用 u′v′Tu'v'^T 算出来的答案就是倍率，求 vv 的过程就是钦定一个对秩没有贡献的原矩阵行对应的元作为自由元，然后高斯消元即可。
还有一些边界情况需要考虑，比如那个“对秩没有贡献的原矩阵行”是全 00 的情况。这种情况的话那么包含这一行的子式肯定全都是 00 了，在我的代码中会消出全 00 的答案，所以正好可以求对。
感觉我写得常数巨大，但还是获得了目前的最优解。以下代码经过了对拍，感觉应该没什么问题了：
CPP
#include<bits/stdc++.h>
#define int long long
using namespace std;
typedef long double ld;
const int N=40,mod=1e9+7,i2=5e8+4;
const ld eps=1e-12;
int T,n,a[N][N],b[N][N],w[N],ans,cnt;
ld A[N][N],c[N][N];
inline ld det(int n){
	ld res=1;
	for(int i=1;i<=n;i++){
		int u=i;
		for(int j=i+1;j<=n;j++)if(fabsl(A[j][i])>fabsl(A[u][i]))u=j;
		if(fabsl(A[u][i])<eps)return 0;
		if(u!=i)res*=-1,swap(A[i],A[u]);
		res*=A[i][i];
		for(int j=i+1;j<=n;j++){
			const ld t=A[j][i]/A[i][i];
			for(int k=i;k<=n;k++)A[j][k]-=A[i][k]*t;
		}
	}
	return res;
}
signed main(){
	ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
	cin>>T;
	while(T--){
		cin>>n,cnt=0,memset(w,0,sizeof(w));
		for(int i=0;i<=n;i++)for(int j=1;j<n;j++)cin>>b[i][j];
		memcpy(a,b,sizeof(b));
		for(int i=1;i<=n;i++)for(int j=n-1;j;j--)a[i][j]^=a[0][j];
		for(int i=1;i<n;i++)for(int j=1;j<n;j++)b[i][j]^=b[n][j],A[i][j]=b[i][j];
		ans=(int)roundl(fabsl(det(n-1)))%mod;
		for(int i=1;i<=n;i++)for(int j=n-1;j;j--)A[j][i]=a[i][j];
		int useless=0;
		for(int i=1,t=1;i<n&&t<=n;i++,t++){
			int u;
			while(t<=n){
				u=i;
				for(int j=i+1;j<n;j++)if(fabsl(A[j][t])>fabsl(A[u][t]))u=j;
				if(fabsl(A[u][t])<eps)useless=t,t++;
				else{cnt++;break;}
			}
			if(t>n)break;
			if(u!=i)swap(A[i],A[u]);
			for(int j=i+1;j<n;j++){
				const ld d=A[j][t]/A[i][t];
				for(int k=t;k<=n;k++)A[j][k]-=A[i][k]*d;
			}
		}
		if(!useless)useless=n; 
		if(cnt==n-1){
			for(int j=1;j<n;j++){
				for(int i=1,t=0;i<=n;i++)if(i!=useless)c[j][++t]=a[i][j];
				c[j][n]=a[useless][j];
			}
			for(int i=1;i<n;i++){
				int u=i;
				for(int j=i+1;j<n;j++)if(fabsl(c[j][i])>fabsl(c[u][i]))u=j;
				if(u!=i)swap(c[i],c[u]);
				for(int j=i+1;j<=n;j++)c[i][j]/=c[i][i];
				c[i][i]=1;
				for(int j=1;j<n;j++){
					if(j==i)continue;
					const ld t=c[j][i];
					for(int k=i;k<=n;k++)c[j][k]-=c[i][k]*t;
				}
			}
			ld v[N];
			v[useless]=-1;
			for(int i=1,t=0;i<=n;i++)if(i!=useless)t++,v[i]=c[t][n];
			for(int i=1,t=0;i<=n;i++)if(i!=useless){
				t++;
				for(int j=1;j<n;j++)A[t][j]=a[i][j];
			} 
			ld t=roundl(det(n-1));
			for(int i=1;i<=n;i++)(ans+=llabs((int)roundl(fabsl(v[i]*t))))%=mod;
		}
		(ans*=i2)%=mod;
		cout<<ans<<"\n";
	}
	return 0;
}

评论

共 0 条评论，欢迎与作者交流。

正在加载评论...