題面

這是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 7

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

/*
          _           _    _            _
    /\   | |         | |  | |          (_)
   /  \  | | _____  _| |__| | ___  _ __ _ _ __   __ _
  / /\ \ | |/ _ \ \/ /  __  |/ _ \| '__| | '_ \ / _` |
 / ____ \| |  __/>  <| |  | | (_) | |  | | | | | (_| |
/_/    \_\_|\___/_/\_\_|  |_|\___/|_|  |_|_| |_|\__, |
                                                 __/ |
                                                |___/
*/