题解：P7429 [THUPC 2017] 气氛

· · 题解

提供一个 O(Tn^3) 的爆标做法。

首先结论是 n−1 维空间中 n+1 个点构成的凸包的体积等于所有选 n 个点构成的凸包体积之和的一半，其他题解对此已有详细说明。

然后考虑 n−1 维空间中 n 个点构成的凸包体积怎么算。先随便选出一个点，然后用其他 n-1 个点到这个点的向量求行列式的绝对值即为体积。

首先根据 OEIS A003432，35 阶所有元素全为 0,1 的行列式的上界是 2\times 3^{36}，证明我还不会。这个说明我们使用 C++ 浮点数类型在本题数据范围下是不会爆精度的。

我们可以把求体积过程中求向量的“减”改为“异或”，这样就是全 0,1 行列式，而且也是对的。考虑这就是把空间转了一下然后把那个点转到原点，而我们只需要知道有向体积的绝对值，所以答案不变。

我们先 O(n^3) 直接计算第 2,3,...n+1 个点组成的凸包体积，然后把除了第一个向量以外的向量全部“异或”上第一个向量。

考虑下面这组数据：

1
4
1 0 1 
0 0 1 
1 1 1 
0 1 1 
0 1 0

“异或”上第一个向量后得到的四个向量是：

1 0 0
0 1 0
1 1 0
1 1 1

我们要求的就是分别去掉每一行后的行列式之和。这个形式有点像代数余子式，我们考虑给这个矩阵左边添上一列：

0 1 0 0
0 0 1 0
0 1 1 0
0 1 1 1

求一个方阵的所有代数余子式是可以做到 O(n^3) 的（首先你需要了解伴随矩阵的概念，这里记为 \operatorname{adj}(A)），过程如下：

接下来就是 \operatorname{rank}(A)=n-1 的情况了。

首先当 \operatorname{rank}(A)=n-1 时 \operatorname{rank}(\operatorname{adj}(A))=1，证明考虑，首先显然 \operatorname{rank}(\operatorname{adj}(A))\ge1，并且有 A\operatorname{adj}(A)=0。然后有矩阵秩经典性质：对于任意矩阵 A_{m\times n} 和 B_{n\times s}，AB=0 时 \operatorname{rank}(A)+\operatorname{rank}(B)\le n，所以 \operatorname{rank}(\operatorname{adj}(A))\le1。

所以 \operatorname{adj}(A) 可以被表示为 uv^T 的形式，其中 u,v 是 n 阶列向量。

因为 A\operatorname{adj}(A)=\operatorname{adj}(A)A=0，所以 u,v 分别在 A 的零空间、左零空间当中，即 Au=A^Tv=0。我们只需要从对应的解空间中分别解出任意一组非零解 u',v'，那么容易发现 uv 是 u'v' 的常数倍。

在本题中已经在左边添上了全零的一列，所以可以直接取 u'=(1,0,0,...,0)。由于本题中我们最后只需要知道绝对值，所以可以省去很多关于 -1 的若干次方的讨论，那么容易发现我们最终只需要求 v'，然后取一个 v' 中的非零位置用朴素方法求一下那个位置对应的子式再除以用 u'v'^T 算出来的答案就是倍率，求 v 的过程就是钦定一个对秩没有贡献的原矩阵行对应的元作为自由元，然后高斯消元即可。

还有一些边界情况需要考虑，比如那个“对秩没有贡献的原矩阵行”是全 0 的情况。这种情况的话那么包含这一行的子式肯定全都是 0 了，在我的代码中会消出全 0 的答案，所以正好可以求对。

感觉我写得常数巨大，但还是获得了目前的最优解。以下代码经过了对拍，感觉应该没什么问题了：

#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;
}