标签：mathbf Point 题解 矩阵 Transform 线性代数 setp double ep
题面
这是2021年CCPC东北四省赛的J题。这里放上CodeForces的链接：J. Transform。下面搬运一下原题面：
Given you two points \((A,B,C)\) and \((x, y, z)\) in 3D space. Let \(L\) be the line cross the point \((A, B, C)\) and the original point \((0, 0, 0)\). You will get two points \(P\) and \(Q\) if you rotate point \((x, y, z)\) around line \(L\) by \(r\) degree and \(-r\) degree. If the \(z\) coordinate of \(P\) is greater than \(Q\), output \(P\). Otherwise, output \(Q\). We guarantee that the solution is unique.
input
This problem contains multiple test cases.
The first line of the input contains an integer \(T (1 \leq T \leq 50000)\), indicating the number of test cases.
Each of the following \(T\) lines contains seven integers \(A, B, C, x, y, z, r\).
\((1 \leq A, B, C, x, y, z \leq 100, 1 \leq r < 180)\).
output
For each test case, output one line contains three real numbers indicates the coordinate of the answer.
Your answer will be accepted if the absolute or relative error between your answer and the standard answer is less than \(10^{-6}\).
sample input
1 1 2 3 4 5 6 7sample output
4.084934830 4.801379781 6.104101869
大意
给出一条通过原点与点\((A,B,C)\)的直线，求给定点\(p=(x,y,z)\)绕该直线顺时针与逆时针旋转\(r\)°（角度制）得到的两个点中\(z\)坐标最大的点的坐标，精确到小数点后6位。
题解
三维空间的旋转问题很难用一个式子来表示，当然可能有现成的结论，但是本人才疏学浅不知道，希望了解的大佬在评论区互相交流~
自然地，我希望能把其转化为二维平面上的旋转问题，这样便可以利用复数的几何意义来求解。如果给定的直线与某条坐标轴重合，例如x轴，那么待旋转的点就相当于在\(yOz\)平面内进行旋转。所以对于任意给出的一条直线，是否有办法将其转化为上述讨论的情况呢？答案是可以的：使用坐标变换将给定的直线转化为x轴，在这个新坐标系下进行旋转，然后再转换回原坐标系。坐标的转换需要什么工具？答案是：线性代数。
我们设给定的直线确定一根向量\(\mathbf{a_1}=(A,B,C)^T\)，很容易找到两条不共线的向量\(\mathbf{a_2, a_3}\)使得这三条向量不共面（见后文）。然后，我们对这三个向量进行施密特正交化，并且进行单位化。定义两个向量\(\mathbf{a,b}\)的内积为\((\mathbf{a,b})\)，则三维向量的施密特正交化过程如下所示：
\[\left\{ \begin{array}{lr} \mathbf{b_1}=\mathbf{a_1},\\ \mathbf{b_2}=\mathbf{a_2}-\frac{(\mathbf{b_1,a_2})}{(\mathbf{b_1,b_1})} \mathbf{b_1},\\ \mathbf{b_3}=\mathbf{a_3}-\frac{(\mathbf{b_1,a_3})}{(\mathbf{b_1,b_1})} \mathbf{b_1}-\frac{(\mathbf{b_2,a_3})}{(\mathbf{b_2,b_2})} \mathbf{b_2} \end{array} \right. \]正交化过程以\(\mathbf{a_1}\)为基准，就是说这个向量的方向是不会变的，这就保证了在线性变换以后它会被变换为坐标轴。然后，我们将正交化的三个列向量单位化得到\(\mathbf{e_1,e_2,e_3}\)，以后按顺序拼接在一起，就得到了由当前坐标系变换到基\(\mathbf{e_1,e_2,e_3}\)张成的空间的过渡矩阵\(C=\mathbf{[e_1,e_2,e_3]}\)。这个过渡矩阵是正交矩阵，可逆，并且在变换前后的线性空间不会发生缩放。求出逆矩阵为\(C^{-1}\)（见后文），那么变换后的给定点坐标\(p'=(x',y',z')\)满足如下式子：
\[(x',y',z')^T=C^{-1}(x,y,z)^T \\ (x,y,z)^T=C(x',y',z')^T \]变换以后就把原问题转换为了点绕x轴旋转的问题，并且由于正交变换不改变形状与相对位置的性质，变换前后旋转的角度是不变的。
在平面内的旋转问题可以利用复数乘法的几何意义：两个复数相乘，辐角相加、模长相乘。记复数\(c=a+b\verb|i| = (a,b)\)，那么定义复数乘法为\((a,b)\times (c,d)=(ac-bd,ad+bc)\)，构造复数\(c_1=(\cos(r),\sin(r)),c_2=(\cos(r),-\sin(r)),c_t=(y',z')\)，分别计算\(c_{t1}=c_1 \times c_t=(y'_1,z'_1),c_{t2} = c_2\times c_t=(y'_2,z'_2)\)，那么可以得到两个点\(p'_1=(x',y'_1,z'_1),p'_2=(x',y'_2,z'_2)\)。这两个点就是基\(\mathbf{e_1,e_2,e_3}\)下的点\(p'\)绕x轴顺逆时针旋转\(r\)°得到的点。接着通过过渡矩阵转换为基\((1,0,0)(0,1,0)(0,0,1)\)下的坐标，就是题面要求的两个点\(P,Q\)。接下来只要判哪个点的竖坐标更大，将其输出即可。
接下来我们解决一些在叙述过程中未解决的问题。
我们知道，在三维坐标系下，如果一个向量的z坐标不为0，那么它必定与\(xOy\)平面不共面，那么要求三条不共面向量，只需要取\(\mathbf{a_2}=(1,0,0),\mathbf{a_3}=(0,1,0)\)即可。另外两种情况同理，这里不再赘述。
要求一个矩阵的逆矩阵，考虑到三维空间的矩阵大小为\(3\times 3\)，在代码实现中可以使用初等变换法，构造矩阵\([C,E]\)，其中\(E\)为单位矩阵，然后对这个矩阵进行初等行变换直到左部变成单位矩阵，即变换为\([E,C^{-1}]\)，这时右部就是我们需要的逆矩阵。
本题给出的时间限制为4000ms，时间非常宽裕（大模拟类题目与计算几何类题目不会在时间上卡你的程序，能做出来的方法一般都会过），并且本题的维数固定为3，每组数据的规模可以认为是一样的，因此不再进行时间复杂度的分析。实际上对于维数\(n\)的复杂度为\(O(n^3)\)，本题\(n=3\)，差别仅在于常数。
下面是代码实现：
#include <bits/stdc++.h>
#define GRP int T;cin>>T;rep(C,1,T)
#define FAST ios::sync_with_stdio(false);cin.tie(0);
#define VOID
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define rrep(i,a,b) for(int i=a;i>=b;--i)
#define elif else if
#define mem(arr,val) memset(arr,val,sizeof(arr))
typedef long long ll;
typedef unsigned long long ull;
const double PI = acos(-1);
inline bool is0(double &q) {
return abs(q) < 1e-10;
} //实数判0：绝对值小于一个极小值
struct Complex { //将复数封装，方便使用
double a, b;
Complex() = default;
Complex(double a, double b): a(a), b(b) {}
void setc(double aa, double bb) {
a = aa, b = bb;
}
Complex operator*(Complex &e) {
return Complex(a * e.a - b * e.b, a * e.b + b * e.a);
} //定义复数乘法
};
struct Point { //封装点（向量）类
double x, y, z;
Point() = default;
Point(double x, double y, double z): x(x), y(y), z(z) {}
double dotc(Point &a) { //点乘
return x * a.x + y * a.y + z * a.z;
}
Point operator*(double a) { //数乘
return Point(x * a, y * a, z * a);
}
Point operator+(Point a) { //加法
return Point(x + a.x, y + a.y, z + a.z);
}
Point operator-(Point a) { //减法
return *this + a * (-1.0);
}
void setp(double x, double y, double z) {
this->x = x;
this->y = y;
this->z = z;
}
void unit() { //单位化
double l = x * x + y * y + z * z;
l = sqrt(l);
x /= l;
y /= l;
z /= l;
}
Point tran(double m[3][3]) { //矩阵乘法
double g[3];
rep(i, 0, 2) {
g[i] = m[i][0] * x + m[i][1] * y + m[i][2] * z;
}
return Point(g[0], g[1], g[2]);
}
};
ostream& operator<<(ostream& os, Point p) { //输出
os << fixed << setprecision(7) << p.x << ' ' << p.y << ' ' << p.z ;
return os;
}
typedef Point Vector;
double a, b, c, x, y, z, deg;
Vector t[3], ep[3];
Point p;
Point pt1, pt2;
double cmx[3][3], cmxr[3][3];
void rev() { //求逆矩阵
double mt[3][6];
double tmp[6];
//拼接矩阵
rep(i, 0, 2) {
rep(j, 0, 2) {
mt[i][j] = cmx[i][j];
}
}
rep(i, 0, 2) {
rep(j, 3, 5) {
if (j - 3 == i) {
mt[i][j] = 1;
} else {
mt[i][j] = 0;
}
}
}
double k;
//初等变换
rep(l, 0, 2) {
if (is0(mt[l][l])) { //对角线元素必须非0
rep(i, l + 1, 2) {
if (!is0(mt[i][l])) {
rep(j, 0, 5) {
tmp[j] = mt[l][j];
}
rep(j, 0, 5) {
mt[l][j] = mt[i][j];
}
rep(j, 0, 5) {
mt[i][j] = tmp[j];
}
break;
}
}
}
k = 1 / mt[l][l];
rep(i, 0, 5) {
mt[l][i] *= k; //元素化为1
}
rep(i, 0, 2) {
if (i == l) continue;
k = -mt[i][l] / mt[l][l];
rep(j, 0, 5) {
mt[i][j] += (k * mt[l][j]);
}
} //该列其余元素化为0
}
rep(i, 0, 2) {
rep(j, 3, 5) {
cmxr[i][j - 3] = mt[i][j];
}
} //获得逆矩阵
}
int main() {
FAST
GRP {
cin >> a >> b >> c >> x >> y >> z >> deg;
deg = deg * PI / 180; //化为弧度制
t[0].setp(a, b, c);
p.setp(x, y, z);
//构造互不共面的向量组
if (a != 0) {
t[1].setp(0, 0, 1);
t[2].setp(0, 1, 0);
}
elif(b != 0) {
t[1].setp(0, 0, 1);
t[2].setp(1, 0, 0);
} else {
t[1].setp(1, 0, 0);
t[2].setp(0, 1, 0);
}
//正交化
ep[0] = t[0];
ep[1] = t[1] - ep[0] * (ep[0].dotc(t[1]) / ep[0].dotc(ep[0]));
ep[2] = t[2] - ep[0] * (ep[0].dotc(t[2]) / ep[0].dotc(ep[0])) - ep[1] * (ep[1].dotc(t[2]) / ep[1].dotc(ep[1]));
//单位化并构造过渡矩阵
rep(i, 0, 2) {
ep[i].unit();
cmx[0][i] = ep[i].x;
cmx[1][i] = ep[i].y;
cmx[2][i] = ep[i].z;
}
rev(); //求逆矩阵
p = p.tran(cmxr); //点的坐标转换
//构造复数并计算
Complex c1(cos(deg), sin(deg)), c2(cos(deg), -sin(deg)), ct(p.y, p.z);
c1 = c1 * ct, c2 = c2 * ct;
pt1.setp(p.x, c1.a, c1.b);
pt2.setp(p.x, c2.a, c2.b);
pt1 = pt1.tran(cmx); //转换回原坐标
pt2 = pt2.tran(cmx);
if (pt1.z > pt2.z) { //输出结果
cout << pt1 << endl;
} else {
cout << pt2 << endl;
}
}
return 0;
}
/*
_ _ _ _
/\ | | | | | | (_)
/ \ | | _____ _| |__| | ___ _ __ _ _ __ __ _
/ /\ \ | |/ _ \ \/ / __ |/ _ \| '__| | '_ \ / _` |
/ ____ \| | __/> <| | | | (_) | | | | | | | (_| |
/_/ \_\_|\___/_/\_\_| |_|\___/|_| |_|_| |_|\__, |
__/ |
|___/
*/
标签：mathbf,Point,题解,矩阵,Transform,线性代数,setp,double,ep 来源： https://www.cnblogs.com/AlexHoring/p/16293080.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。