2 条题解
-
0
题目大意
你在一个 维空间中。
给出空间中 个点的坐标 及它们与你的距离 ，
请你求出你所在的坐标 。
若有多个解，输出任意一个。
，坐标范围 ，保证数据随机且有解。
题解
给出一个 naive 的数形结合的做法。
根据题意，不难列出方程
$$\begin{cases} (\vec{x}-\vec{y}_1)^2=r_1^2\\ (\vec{x}-\vec{y}_2)^2=r_2^2\\ \cdots \\ (\vec{x}-\vec{y}_n)^2=r_n^2 \end{cases}$$这个方程含有平方，所以没法直接高斯消元。
但注意到
$$(\vec{x}-\vec{y})^2=\vec{x}^2-2\vec{y}\cdot\vec{x}+\vec{y}^2$$其中 是相同的。
考虑将原点平移到 ，设 $\vec{x}'=\vec{x}-\vec{y}_1,\forall k,\vec{y}_k'=\vec{y}_k-\vec{y}_1$，则
$$\begin{cases} \vec{x}'^2=r_1^2&(1)\\ (\vec{x}'-\vec{y}_2')^2=r_2^2\\ (\vec{x}'-\vec{y}_3')^2=r_3^2\\ \cdots\\ (\vec{x}'-\vec{y}_n')^2=r_n^2 \end{cases}$$下面 个式子同时减去 ，得到
$$\begin{cases} -2\vec{y}'_2\cdot\vec{x}+\vec{y}_2'^2=r_2^2-r_1^2\\ -2\vec{y}'_3\cdot\vec{x}+\vec{y}_3'^2=r_3^2-r_1^2\\ \cdots\\ -2\vec{y}'_n\cdot\vec{x}+\vec{y}_n^2=r_n^2-r_1^2 \end{cases}$$这等价于线性方程组
$$A=\begin{bmatrix} (\vec{y}_2')^T\\ \hline(\vec{y}_3')^{T^{^{^{}}}}\\ \hline\vdots \\ \hline(\vec{y}_n')^{T^{^{^{}}}} \end{bmatrix},\vec{a}=-\dfrac{1}{2}\begin{bmatrix} r_2^2-r_1^2-\vec{y}_2^2\\ r_3^2-r_1^2-\vec{y}_3^2\\ \vdots\\ r_n^2-r_1^2-\vec{y}_n^2 \end{bmatrix}$$高斯消元求解该线性方程组，若其有唯一解，则解就是答案。
如果解不是唯一的怎么办？这时“形”便派上用场。
让我们转换一下题意：给定 维空间中 个超球，
球心为 ，半径为 ，求它们的交点。
先把这个问题搁置一边，来看看这个高中文化课常用的 trick：
给定二维空间两个圆的方程：
$$\begin{cases} (x-a_1)^2+(y-b_1)^2=r_1^2&(1)'\\ (x-a_2)^2+(y-b_2)^2=r_2^2&(2)' \end{cases}$$它们的交线如图所示：
现在我们要求交线的方程，一般人的想法是将 联立，求得交点后再求交线。
但最快的方法是将 相减，直接求得交线的方程
$$(x-a_1)^2-(x-a_2)^2+(y-b_1)^2-(y-b_2)^2=r_1^2-r_2^2$$接下来我们简短证明下这么做的正确性：
- 中 系数均为 ，因此 一定是个二元一次不定方程，
也就一定是一条直线。 - ，因此圆 与该直线的交点一定在圆 上。
同理，圆 与该直线的交点一定在圆 上。
因此该直线与两个圆中一个的交点等价于两个圆的交点。
综上， 就是两个圆的交线方程。
显然这个结论及其证明可以推广到高维空间，若已知三维空间两个球的方程：
$$\begin{cases} (x-a_1)^2+(y-b_1)^2+(z-c_1)^2=r_1^2&(1)''\\ (x-a_2)^2+(y-b_2)^2+(z-c_2)^2=r_2^2&(2)''\\ \end{cases}$$它们的交面如图所示：
把 和 相减，可直接得到交面的方程
$$(x-a_1)^2-(x-a_2)^2+(y-b_1)^2-(y-b_2)^2+(z-c_1)^2+(z-c_2)^2=r_1^2-r_2^2$$且交面与其中一个球的交点等价于两球的交点。
同理，若已知 维空间的两个超球：
$$\begin{cases} (\vec{x}-\vec{y}_1)^2=r_1^2\\ (\vec{x}-\vec{y}_2)^2=r_2^2 \end{cases}$$将两个超球的方程相减，可得到 维超交面的方程
$$(\vec{x}-\vec{y}_1)^2-(\vec{x}-\vec{y}_2)^2=r_1^2-r_2^2$$且超交面与其中一个超球的交点等价于两个超球的交点。
此时我们已经能解释之前的高斯消元在做什么了：
- 将所有方程减去 等价于求超球 与其他所有超球的超交面方程。
- 高斯消元等价于求这些超交面的交。
- 这些超交面的交与超球 的交点都是合法的解，即所有超球的交点。
分析完这些后，解不唯一的情况变成了一个这样的问题：
设高斯消元后自由元的个数为 ，
则我们需要求解一个 维超球和一个 维 “平面” （准确的说是线性流形）的交点。
同样考虑二维空间的情况：求解一个圆和一条直线的交点。
一个显然的做法是过圆心向直线作垂线，然后勾股定理算出垂足到交点的距离：
拓展到三维空间时，情况是类似的：
-
对于三维空间一条直线和球相交，跟刚才的情况是一样的。
-
对于三维空间一个平面和球相交，可以确定截面一定是个圆。 同样作垂线得到圆心后勾股定理求解。
想一想怎么归纳到高维空间，提示一下：一维空间的超球是两个点。
答案揭晓： 维空间一个超球与一个 维线性流形相交，截面一定是个 维超球。
过原超球球心向线性流形作垂线即可得到截面的球心，再勾股定理即可得到球心的半径。我们由此得到求解原问题的一个大致的思路：
- 过超球心向线性流形作垂线，求出垂足的坐标 。
- 勾股定理算出截面的半径 。
- 在线性流形上随便找一个长度为 的向量 ， 就是一个合法的解。
回到之前的线性方程组。高斯消元后，其一定会变为如下形式：
$$\left[ \begin{array}{ccccc|} 1&&&&&v_{1,1}&v_{2,1}&\cdots&v_{k,1}\\ &1&&&&v_{1,2}&v_{2,2}&\cdots&v_{k,2}\\ &&\ddots&&&\vdots&\vdots&\vdots&\vdots\\ &&&1&&v_{1,d-k-1}&v_{2,d-k-1}&\cdots&v_{k,d-k-1}\\ &&&&1&v_{1,d-k}&v_{2,d-k}&\cdots&v_{k,d-k}\\ \hline 0&0&0&0&0&0&0&0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots \end{array}\right] \begin{bmatrix} x_1\\x_2\\\vdots\\x_{d-k}\\\hline x_{d-k+1}\\ x_{d-k+2}\\ \vdots \\x_d \end{bmatrix}=\begin{bmatrix} a_1'\\ a_2'\\ \vdots\\ a_{d-k}'\\\hline 0\\0\\ \vdots \\0 \end{bmatrix}$$简记为如下形式：
$$\left[\begin{array}{c|c|c|c|c} I&\vec{v}_1&\vec{v}_2&\cdots&\vec{v}_k\\ \hline O&\vec{0}&\vec{0}&\vec{0}&\vec{0} \end{array}\right] \begin{bmatrix}\vec{x}_0\\ \hline\vec{x}_1\end{bmatrix} =\begin{bmatrix}\vec{a}'\\ \hline\vec{0}\end{bmatrix}$$我们需要确定线性方程组的解集，从而确定整个线性流形。
- 若 ，显然 。
因此 $\begin{bmatrix}\vec{a}'\\ \hline\vec{0}\end{bmatrix}$ 是该线性方程组的一组特解。 - 若 ，
则 $\vec{x}_0+\lambda\vec{v}_1=\vec{a}',\vec{x}_0=\vec{a}-\lambda\vec{v}_1$。
因此 每增加 ， 就要减去 。
综上，设
$$\vec{u}_1=\begin{bmatrix}-\vec{v}_1\\\hline 1\\0\\0\\0\\ \vdots\end{bmatrix} ,\vec{u}_2=\begin{bmatrix}-\vec{v}_2\\\hline 0\\1\\0\\0\\ \vdots\end{bmatrix}, \vec{u}_3=\begin{bmatrix}-\vec{v}_3\\\hline 0\\0\\1\\0\\ \vdots\end{bmatrix},\cdots$$$$\vec{b}=\begin{bmatrix} \vec{a}'\\ \hline \vec{0} \end{bmatrix},V=\{[\vec{u}_1|\vec{u}_2|\cdots|\vec{u}_k]\vec{x}_1|\vec{x}_1\in\mathbb{R}_k\}$$则线性方程组的解集为
即线性空间 内的所有向量加上特解 ，与线性流形的定义一致。
接下来考虑作垂线，由于球心为 ， 的基底为 ，
$$\vec{u}_1\cdot \vec{p}=\vec{u}_2\cdot\vec{p}=\cdots=\vec{u}_k\cdot \vec{p}=0 \qquad{(2)}$$
我们的目标是找出 ，使得 与 中所有的基向量垂直，即把 拆开：
$$\vec{p}=\vec{b}+[\vec{u}_1|\vec{u}_2|\cdots|\vec{u}_k]\vec{x}_1\qquad{(3)}$$然后将 表示成矩阵形式：
$$\begin{bmatrix} \vec{u}_1^T\\ \hline\vec{u}_2^{T^{^{^{}}}}\\ \hline\cdots\\ \hline\vec{u}_k^{T^{^{^{}}}} \end{bmatrix}(\vec{b}+[\vec{u}_1|\vec{u}_2|\cdots|\vec{u}_k]\vec{x}_1)=\vec{0}$$注意到 $\begin{bmatrix}\vec{u}_1^T\\\hline\vec{u}_2^{T^{^{^{}}}}\\\hline\cdots\\\hline\vec{u}_k^{T^{^{^{}}}} \end{bmatrix}[\vec{u}_1|\vec{u}_2|\cdots|\vec{u}_k]$ 为 阶方阵 ，因此高斯消元解线性方程组
$$\begin{bmatrix}\vec{u}_1^T\\\hline\vec{u}_2^{T^{^{^{}}}}\\\hline\cdots\\\hline\vec{u}_k^{T^{^{^{}}}} \end{bmatrix}[\vec{u}_1|\vec{u}_2|\cdots|\vec{u}_k]\vec{x}_1=-\begin{bmatrix}\vec{u}_1^T\\\hline\vec{u}_2^{T^{^{^{}}}}\\\hline\cdots\\\hline\vec{u}_k^{T^{^{^{}}}} \end{bmatrix}\vec{b}$$即可得到 ，代入 即可求出 。
因为过线性流形（直线/平面）上一点有且只有一条直线与该流形（直线/平面）垂直，
这个线性方程组一定有唯一解，也就一定能解出 。
解出 后，勾股定理得到
然后 就是一组合法的解了。
最后不要忘记，我们一开始将原点平移到了 ，要平移回去才行哦！
Code
/* this code is made by warzone 2022.3.11 15:30 */ #include<stdio.h> #include<string.h> #include<math.h> #include<algorithm> typedef double db; struct READ{//快读 char c,w; inline READ(){c=getchar();} template<typename type> inline READ& operator >>(type& num){ for(w=1;'0'>c||c>'9';c=getchar()) w=c=='-'? -1:1; for(num=0;'0'<=c&&c<='9';c=getchar()) num=num*10+(c-'0'); if(c=='.') for(db i=1;c=getchar(),'0'<=c&&c<='9';) i*=0.1L,num+=i*(c-'0'); return num*=w,*this; } }cin; int d,n; int pos[512],pos_[512]; db y[512],r; db A[512][512],a[512]; db B[512][512],b[512]; int top,top_; template<typename type1,typename type2,typename type3> inline void gauss(type1 &A,type2 &a,type3 &pos,//高斯消元 int &top,const int d,const int n){ top=d; for(int id=0;id<top&&id<n;++id){ while(A[id][id]==0){ for(int i=id+1;i<n;++i) if(A[i][id]){ for(int j=id;j<d;++j) std::swap(A[id][j],A[i][j]); std::swap(a[id],a[i]);break; } if(A[id][id]==0){//把消元失败的维度放到最后面 if(id>=--top) return; std::swap(pos[top],pos[id]); for(int i=0;i<n;++i) std::swap(A[i][id],A[i][top]); } } if(A[id][id]!=1){ db get=A[id][id]; for(int i=id+1;i<d;++i) A[id][i]/=get; A[id][id]=1,a[id]/=get; } for(int i=0;i<n;++i) if(i!=id&&A[i][id]){ db get=A[i][id]; for(int j=id+1;j<d;++j) A[i][j]-=A[id][j]*get; A[i][id]=0,a[i]-=a[id]*get; } } } inline void print(){//平移回来并输出答案 for(int i=0;i<d;++i) y[pos[i]]+=a[i]; for(int i=0;i<d;++i) printf("%0.6lf ",y[i]); } int main(){ cin>>d>>n; for(int i=0;i<d;++i) cin>>y[i],pos[i]=i;//输入 (1) 式 cin>>r,--n; for(int id=0;id<n;++id){ for(int i=0;i<d;++i){//减去 (1) 式 cin>>A[id][i],A[id][i]-=y[i]; a[id]+=A[id][i]*A[id][i]; } db get; cin>>get,a[id]=(a[id]+r*r-get*get)/2; } gauss(A,a,pos,top,d,n),top=std::min(top,n);//第一次高斯消元 if(top==d) return print(),0;//唯一解 for(int i=top;i<d;++i) A[i][i]=-1; const int k=d-top; for(int id=0;id<k;++id){//计算 \vec{p} 的线性方程组 pos_[id]=id; for(int i=0;i<k;++i) for(int j=0;j<d;++j) B[id][i]+=A[j][id+top]*A[j][i+top]; for(int j=0;j<d;++j) b[id]-=A[j][id+top]*a[j]; } gauss(B,b,pos_,top_,k,k);//第二次高斯消元 for(int id=0;id<k;++id)//将 \vec{x}_1 代回 \vec{p} for(int i=0;i<d;++i) a[i]+=b[id]*A[i][id+top]; db r_=r*r,r__=0; for(int i=0;i<d;++i){//勾股定理计算 r' r__+=A[i][top]*A[i][top]; r_-=a[i]*a[i]; } r_=sqrtl(r_/r__); for(int i=0;i<d;++i) a[i]+=A[i][top]*r_; return print(),0; } - 中 系数均为 ，因此 一定是个二元一次不定方程，
-
0
备注：
\pmb在题解区可以正常渲染。仅预览不支持。Update on 2.24：规范数学公式用法（去掉了所有矩阵前面的
\pmb）。Update on 4.14：修改一个 Typo。
P8141 [ICPC2020 WF] What’s Our Vector, Victor?
给定 个 维向量 与 ，求 维向量 满足 。保证有解。。
根据题意，可以列出以下方程。
$$\begin{cases} (x_1 - v_{1, 1}) ^ 2 + (x_2 - v_{1, 2}) ^ 2 + \cdots = e_1 ^ 2 \\ (x_1 - v_{2, 1}) ^ 2 + (x_2 - v_{2, 2}) ^ 2 + \cdots = e_2 ^ 2 \\ \cdots \end{cases}$$将平方项展开后，出现了 项。众所周知高斯消元只能处理线性方程组，所以考虑常用套路 差分消去平方项。处理后，我们得到了 个 元线性方程组，和最开始的一个 元二次方程。形如
$$\begin{cases} (x_1 - v_{1, 1}) ^ 2 + (x_2 - v_{1, 2}) ^ 2 + \cdots = e_1 ^ 2 \\ 2(v_{1, 1} - v_{2, 1}) x_1 + 2(v_{1, 2} - v_{2, 2}) x_2 + \cdots = e_2 ^ 2 - e_1 ^ 2 + v_{1, 1} ^ 2 - v_{2, 1} ^ 2 + v_{1, 2} ^ 2 - v_{2, 1} ^ 2 + \cdots \\ 2(v_{2, 1} - v_{3, 1}) x_1 + 2(v_{2, 2} - v_{3, 2}) x_2 + \cdots = e_3 ^ 2 - e_2 ^ 2 + v_{2, 1} ^ 2 - v_{3, 1} ^ 2 + v_{2, 2} ^ 2 - v_{3, 1} ^ 2 + \cdots \\ \cdots \end{cases}$$用向量和矩阵描述 前的系数，则有
$$\begin{cases} \| \pmb x - \pmb v_1 \| ^ 2 = e_1 ^ 2 \\ A\pmb x = \pmb b \end{cases}$$其中 为 矩阵， 为 维向量即所求答案， 为 维向量。
首先考虑化简二次式。令 。则方程可写为
$$\begin{cases} \| \pmb x' \| ^ 2 = e_1 ^ 2 \\ A\pmb x' = \pmb b - A\pmb {v}_1 \end{cases}$$接下来我们用 代替 ，求解后再加上 。
高斯消元求解线性方程组 ，得到线性方程组的解集的 参数向量形式。即 $\pmb x = \pmb x_{\rm base} + \sum\limits_{i = 1} ^ r c_i\pmb v_i$，其中 是 自由元 个数。
-
众所周知，对于线性方程组的解集，任何 主元 都可以写作自由元的 线性组合 形式。考虑令所有自由元等于 ，我们得到了 。将所有自由元移到方程右侧，我们得到了所有主元关于所有自由元的线性组合的权，稍做变形即可得到对应向量。
例如一个简化行阶梯型增广矩阵对应的线性方程组
$$\begin{cases} x_1 + 6x_2 + 3x_4 = 0 \\ x_3 - 4x_4 = 5 \\ x_5 = 7 \end{cases}$$其中 是主元， 是自由元。因此，其解集的参数表示为
$$\begin{cases} x_1 = -6x_2 - 3x_4 \\ x_2 = x_2 \\ x_3 = 5 + 4x_4 \\ x_4 = x_4 \\ x_5 = 7 \end{cases}$$即 $\pmb x = \begin{bmatrix} -6x_2 - 3x_4 \\ x_2 \\ 5 + 4x_4 \\ x_4 \\ 7 \end{bmatrix}$。将常数， 和 分离，得到 $\pmb x = \begin{bmatrix} 0 \\ 0 \\ 5 \\ 0 \\ 7 \end{bmatrix} + x_2\begin{bmatrix} -6 \\ 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} + x_4\begin{bmatrix} -3 \\ 0 \\ 4 \\ 1 \\ 0 \end{bmatrix}$。即方程组的解可写为一个常向量 加上若干向量的线性组合。线性组合的权对应了自由变量的取值。
为了让 满足 ，通常采用的手段是先求出 的 最小范数解 ，再通过调整法求解。显然，如果 则无解，与题意矛盾，因此 。而通过加上足够大的 （任意一个均可），我们总能让得到的 的模长趋于无穷大。又因为 即 的变化随着 的权 的变化 连续，因此调整法正确。
- 关于最小范数解 ，它与所有 垂直。若不与某个 垂直，令 且 ，通过点积的性质和数值分析容易证明 。因此可以列出方程 ，即 。若将 看成矩阵 $V = \begin{bmatrix} \pmb v_1 & \pmb v_2 & \cdots & \pmb v_r \end{bmatrix}$，则 。又因为 $\pmb x = \pmb x_{\rm base} + \sum\limits_{i = 1} ^ r c_i \pmb v_i = \pmb x_{\rm base} + V\pmb c$，其中 是 的列向量，因此 。据此可高斯消元求得 ，继而求出对应的 。
由于 垂直于所有 ，所以调整时直接使用勾股定理，随便选出一个向量 ，然后令 $\pmb x = \pmb x_{\min} + \sqrt {e_1 ^ 2 - \| \pmb x_{\min} \| ^ 2} \times \dfrac {\pmb v_1} {\| \pmb v_1\|}$ 即可。时间复杂度 。
#include <bits/stdc++.h> using namespace std; const int N = 500 + 5; const double eps = 1e-9; int n, d, fr, ma, fre[N], mai[N], rev[N]; double v[N][N], A[N][N], dis[N], ans[N]; double b[N], xm[N], co[N][N]; int find(int m, double *A) { for(int i = 1; i <= m; i++) if(fabs(A[i]) > eps) return i; return -1; } template <class T> void Gauss(int n, int m, T A) { int r = 1; for(int c = 1; c < m && r <= n; c++) { int p = r; for(int i = r + 1; i <= n; i++) if(fabs(A[i][c]) > fabs(A[p][c])) p = i; swap(A[r], A[p]); if(fabs(A[r][c]) < eps) continue; for(int i = m; i >= c; i--) A[r][i] /= A[r][c]; for(int i = r + 1; i <= n; i++) for(int j = m; j >= c; j--) A[i][j] -= A[r][j] * A[i][c]; r++; } for(int i = n - 1; i > 0; i--) { for(int j = i + 1; j <= n; j++) { int main = find(m - 1, A[j]); if(main == -1) continue; if(fabs(A[i][main]) > eps) for(int k = m; k >= main; k--) A[i][k] -= A[i][main] * A[j][k]; } int main = find(m - 1, A[i]); if(main == -1) continue; for(int j = m; j >= main; j--) A[i][j] /= A[i][main]; } } void init() { for(int i = 1; i < n; i++) { int main = find(d, A[i]); if(main != -1) rev[main] = -1, mai[++ma] = main; } for(int i = 1; i <= d; i++) if(!rev[i]) fre[++fr] = i, co[rev[i] = fr][i] = 1; } void solve() { for(int i = 1; i < n; i++) { int main = find(d, A[i]); b[main] = A[i][d + 1]; for(int j = main + 1; j <= d; j++) if(~rev[j]) co[rev[j]][main] = -A[i][j]; } memset(A, 0, sizeof(A)); for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) A[i][fr + 1] -= co[i][j] * b[j]; for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) for(int k = 1; k <= fr; k++) A[i][k] += co[i][j] * co[k][j]; Gauss(fr, fr + 1, A); for(int i = 1; i <= fr; i++) for(int j = 1; j <= d; j++) xm[j] += A[i][fr + 1] * co[i][j]; for(int i = 1; i <= d; i++) xm[i] += b[i]; double normx = 0, normc = 0; for(int i = 1; i <= d; i++) normx += xm[i] * xm[i], normc += co[1][i] * co[1][i]; normc = sqrt(normc); for(int i = 1; i <= d; i++) ans[i] = xm[i] + co[1][i] / normc * sqrt(dis[1] * dis[1] - normx); } int main() { cin >> d >> n; for(int i = 1; i <= n; i++) { for(int j = 1; j <= d; j++) scanf("%lf", &v[i][j]); scanf("%lf", &dis[i]); } for(int i = 1; i < n; i++) { A[i][d + 1] = dis[i + 1] * dis[i + 1] - dis[i] * dis[i]; for(int j = 1; j <= d; j++) { A[i][j] = 2 * (v[i][j] - v[i + 1][j]); A[i][d + 1] += v[i][j] * v[i][j] - v[i + 1][j] * v[i + 1][j]; } } for(int i = 1; i < n; i++) for(int j = 1; j <= d; j++) A[i][d + 1] -= A[i][j] * v[1][j]; Gauss(n - 1, d + 1, A), init(); if(fr) solve(); else for(int i = 1; i <= d; i++) ans[i] = A[i][d + 1]; for(int i = 1; i <= d; i++) printf("%.17lf\n", ans[i] + v[1][i]); return 0; } /* stupid mistakes: forget to check if main is -1 注意 co[i][j] 的意义，表示 v_{i_{j, 1}}，即第 i 个自由元改变时，x_j 的变化系数 勾股定理忘记开根 =.= v1 正负搞错 求范数忘记开根 */ -
- 1