BOJ 31383 - Rotation Transformation (Python3)

solved.ac 루비를 달고 나서 처음으로 푸는 문제다. 깡구현으로 풀던지 행렬의 성질을 이용해서 풀던지 어느 쪽이든 테이스트 차이인데, 수학을 한 사람인 내가 보기에는 행렬의 성질을 사용해서 푸는 문제라고 느껴졌다. 재미있게 풀었다. 추천할 법해서 포스트로 풀이를 작성한다.

 

그럼 본격적으로 풀이에 돌입해 보자.

---

문제

A simple motion over and over... - t.A.T.u., <A Simple Motion>

줄리아는 간단한 동작을, 특히 3차원 공간에서의 간단한 동작을 좋아한다. 최근 그녀는 원점을 통과하는 축을 중심으로 회전하는 간단한 동작에 대응하는 행렬을 발견했다. 회전은 회전축을 나타내는 단위벡터 \(v=(v_1, v_2, v_3)^\top\)와 회전 각 \(\alpha\)로 정의된다.

\(v\)의 종점에서 원점을 바라봤을 때, 회전이 반시계방향으로 이루어지면 \(\alpha\)는 양수이고 시계방향으로 이루어지면 음수이다.

회전행렬은 점 \((x,y,z)^\top\)을 \((x^\prime,y^\prime,z^\prime)^\top\)으로, 다음과 같이 움직인다.

\[\begin{pmatrix}x^\prime\\y^\prime\\z^\prime\end{pmatrix}=\begin{pmatrix}A_{11}&A_{12}&A_{13}\\A_{21}&A_{22}&A_{23}\\A_{31}&A_{32}&A_{33}\end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}\]

불행하게도 줄리아는 \(v,\alpha\)가 무슨 값인지 모른다. 그녀는 당신에게 회전행렬로부터 회전축 \(v\)와 회전각 \(\alpha\)를 복원해 달라고 의뢰했다. 복원이 끝나면 그녀는 간단한 동작을 계속 수행할 것이다.

 

입력

회전행렬 \(A\)가 주어진다. 모든 \(i,j\)에 대해서 \(|A^\prime_{ij}A_{ij}|<10^{-13}\)인 회전행렬 \(A^\prime\)이 실제로 존재함이 보장된다.

 

예제 입력)

0 -1 0
1 0 0
0 0 1

 

출력

첫 번째 줄에 60분법으로 나타낸 회전각 \(\alpha\)를, 두 번째 줄에는 회전축이 되는 단위벡터 \(v\)의 각 성분을 출력한다. \(1\le |\alpha|\le 179\)임이 보장된다.

모든 값을 최대한 정확하게 출력하라. 당신이 출력한 회전각과 단위벡터로부터 만들어진 행렬과 주어진 행렬의 모든 성분의 차가 \(10^{-6}\)을 넘지 않는다면 정답으로 인정된다.

 

예제 출력)

90
0 0 1

 

---

내 코드

from math import acos,pi,sin
mat = [[*map(float,input().split())] for _ in range(3)]
trace = sum(mat[i][i] for i in range(3))
angle = acos((trace-1)/2)
sa = sin(angle)
print(angle*180/pi)
print((mat[2][1]-mat[1][2])/(2*sa), (mat[0][2]-mat[2][0])/(2*sa), (mat[1][0]-mat[0][1])/(2*sa))

 

이 풀이 속에서는 \(\alpha\)를 편의상 호도법으로 나타낼 것이다.

회전축과 수직인 두 벡터를 \(x,y\)라고 한다. 그렇다면 \(\{v,w,u\}\)는 orthonormal basis이다.

유클리드 공간의, 좌표축에 평행한 세 단위벡터로 구성된 orthonormal basis \((\hat{i},\hat{j},\hat{k})\)를 \(\{v,w,u\}\)로 보내는 변환을 나타내는 행렬을 \(Q\)라고 하자. 정확하게는 \(Q=(v, w, u)\)이다. 그렇다면 \(Q^{-1}AQ=\begin{pmatrix}1&0&0\\0&\cos\alpha&-\sin\alpha\\0&\sin\alpha&\cos\alpha\end{pmatrix}\)가 된다. \(\alpha\ne n\pi\)이므로 \(Q^{-1}AQ\)는 대칭행렬이 아니므로, \(A\) 역시 대칭행렬이 아니다. (이 사실은 증명하기 쉬우니 직접 한 번 해 보라.)

 

더욱 중요한 사실은, \((Q^{-1}AQ)^\top=(Q^{-1}AQ)^{-1}\)라는 것이다. 이 사실 역시 증명하기 쉬우니 넘어가겠다. 어쨌든, 이 사실을 이용해서 \(\begin{pmatrix}2&0&0\\0&2\cos\alpha&0\\0&0&2\cos\alpha\end{pmatrix}=Q^{-1}(A+A^\top)Q\)를 얻을 수 있다. 여기서 회전각을 얻을 수 있다. 양 변에 행렬의 trace, 대각합을 취하면 다음과 같은 식을 얻는다.

\[2+4\cos\alpha=\text{Tr}\begin{pmatrix}2&0&0\\0&2\cos\alpha&0\\0&0&2\cos\alpha\end{pmatrix}=\text{Tr}(Q^{-1}(A+A^\top)Q)=\text{Tr}(A+A^\top)=2\text{Tr}(A)\]

그러므로 \(A\)의 trace로부터 회전각 \(\alpha=\cos^{-1}\frac{\text{Tr}(A)-1}{2}\)을 복구해 낼 수 있다.

 

이제는 축이 되는 단위벡터를 복구할 타이밍이다. \(Av=v\)는 당연히 성립한다. 회전축에 평행한 벡터는 회전운동에 의해서 변하지 않는다. 따라서 \(A^\top v = A^{-1}v=v\)이 성립한다. 이 두 식을 빼 주면 \((A-A^\top)v=0\)이 성립한다. 따라서, \(v\)는 다음과 같은 행렬의, eigenvalue가 \(0\)인 eigenvector여야 한다는 것이다.

\[\begin{pmatrix}0&A_{12}-A_{21}&A_{13}-A_{31}\\A_{21}-A_{12}&0&A_{23}-A_{32}\\A_{31}-A_{13}&A_{32}-A_{23}&0\end{pmatrix}\]

따라서 \(v\)는 \(e=(A_{23}-A_{32}, A_{13}-A_{31},A_{21}-A_{12})^\top\)에 평행이다. 이 vector의 norm을 직접 구해줘도 되지만, 이 값이 \(2\sin\alpha\)임을 이용하는 방법이 있다. 직접 구해줄 때는, \(v=\pm\|e\|^{-1}e\) 양 쪽 다 가능함을 유의할 필요가 있다. 값을 이용할 때는 정당화가 필요하다면 논문을 참조하자. 정당화가 필요 없으면 적당한 자료조사로 종료할 수 있다. 나는 Efficient conversion from rotating matrix to rotation axis and angle by extending Rodrigues` formula를 참고했다.

 

이로서 31383번의 풀이를 마친다.