단위 쿼터니언의 로그와 지수 함수를 빠르게 계산한다

2026년 10월 10일

지수 공간을 단위 공으로 접었을 때 단위 쿼터니언의 log와 exp를 효율적이고 정확하게 계산하는 특수화 함수를 살펴본다. 목표는 표준 함수를 그대로 쓰는 구현보다 오차 한계가 커지지 않도록 하는 것이다. 예시는 스칼라 연산이지만 SIMD에도 적합하도록 설계했다.

특수화하는 함수는 다음과 같다.

[\begin{align*}

\mathbf{v} & = \frac{2}{\pi}\tfunc{log}{\abs{Q}} \

Q & = \tfunc{exp}{\frac{\pi}{2}~\mathbf{v}}

\end{align*}]

여기서 $\abs{Q}$는 스칼라 성분이 양수가 되도록 정규화한 쿼터니언을 뜻한다. 즉 $w$가 음수이면 $Q$의 부호를 뒤집는다. 이 특수화로 지수 공간은 토크를 최소화하는 각도 범위로 제한된다.

선형 공간에서 지수 공간으로

먼저 표준 라이브러리 함수를 사용해 log를 정확하게 계산하는 방법을 보자. 스케일 인자 $m$과 단위 쌍벡터 $\mathbf{u}$를 갖는 임의의 쿼터니언 $Q$가 다음과 같다고 하자.

[\begin{align*}

Q & = m~\cos\left(\theta\right)+m~\sin\left(\theta\right)~\mathbf{u} \

& = w + \mathbf{v}

\end{align*}]

그러면 $\tfunc{log}{Q} = \theta~\mathbf{u}$를 계산한다. 표준 함수를 사용하는 구현은 다음과 같다.

// log(q) : standard functions static inline vec3f_t quatf_log(quatf_t q) { float x = q[3]; // m cos(Θ) float t = quat_bnorm_fma(q); // m² sin(Θ)² float y = sqrtf(t); // |V| = m sin(Θ) float a = atan2f(y,x); // Θ : this does all the heavy lifting float s = a/y; // scale factor // 'y' can only be "normal" or zero. In the limit the scale factor 's' // approaches one (as y approaches zero). Alternately this could set // s to zero for y==0 which would filter tiny angles to the identity. s = (y != 0.f) ? s : 1.f; return s*quat_bivector(q); }

앞의 함수를 감싸 특수화 버전을 만들 수 있다. 근사 거리 오차는 대략 반 바퀴 각도 오차다.

// wrap log: distance error ~2.54740144e-07 static inline vec3f_t quatf_fem_hq(quatf_t q) { static const float K = 0x1.45f306p-1f; // 2/π q = (q[3] >= 0.f) ? q : -q; return K*quat_log(q); }

여기서 목표는 앞의 구현과 같은 오차 한계를 내는 것이 아니다. 대신 $Q$의 크기가 정확히 1이라고 가정하는 버전을 표준 함수로 만든다. 표준 log 구현은 생략하고 특수화 구현부터 보자.

// (2/π) log(|q|) : standard functions: move to acos and computing |V| as sqrt(1-w²) static inline vec3f_t quatf_fem_std(quatf_t q) { static const float K = 0x1.45f306p-1f; // 2/π float k = copysignf(K,q[3]); // sgn(q.w) 2/π float w = fabsf(q[3]); // |q.w| w = (w <= 1.f) ? w : 1.f; // clamp |w| to range // distance error (approximate) float y = sqrtf(-fmaf(w,w,-1.f)); // 3.40125329e-07 //float y = sqrtf((1.f-w)*(1.f+w)); // 3.63783706e-07 ('y' error bound is same as FMA above) //float y = sqrtf(1.f-w*w); // 5.11455530e-07 //float y = sqrtf(quat_bnorm_fma(q)); // 4.68131871e-06 float a = acosf(w); float s = a/y; s = (y != 0.f) ? s : 1.f; // same comments on limit s *= k; return s*quat_bivector(q); }

코드에는 $y = \abs{\mathbf{v}}$를 계산하는 여러 방법이 주석으로 남아 있다. 여기서는 모든 대상 플랫폼에 하드웨어 FMA가 있다고 가정한다. 이를 지원하지 않는 플랫폼에서는 나머지 경우에 곱셈과 덧셈으로 바꿀 수 있다. 플랫폼마다 비트 단위로 같은 결과를 얻으려면 어느 방식을 쓸지 일관되게 정해야 한다. 그 밖의 경우에는 오차 한계가 조금 커질 뿐이다.

$y$를 계산할 때는 첫 번째 주석 처리된 방법을 쓰는 편이 낫다. 두 번째 방법은 FMA를 곱셈과 덧셈으로 바꾸면 오차가 커지는 것을 보여준다. 마지막 방법은 $y$ 자체를 정확히 구해도 각도 계산에 그 정보가 쓰이지 않아 오히려 가장 나쁜 결과를 낸다. quatf_log에서 float y = sqrtf(-fmaf(w,w,-1.f))를 계산하도록 바꾸면 같은 문제가 생긴다. 감싼 함수 quatf_fem_hq의 오차는 5.11220598e-07까지 커져 위 구현보다 나빠진다.

직접 구현하려면 위 식의 스케일 인자 $s$를 $w$만의 함수로 근사하면 된다.

[\func{s}{w} = \frac{2}{\pi}~\frac{\tfunc{acos}{w}}{\sqrt{1-w^2}}]

근사 방법은 여러 가지다. 파일에는 현재 오차 한계가 서로 다른 19가지 선택지가 있으며, 다항식·유리식·변환 근사의 세 가지 형태를 사용한다. 변환 근사는 정확한 하드웨어 역제곱근이 없으면 큰 의미가 없다. 여기서는 앞서 본 구현보다 오차 한계가 작은 첫 번째 유리 근사를 소개한다.

// ~3.37359953e-07 static inline vec3f_t quatf_fem(quatf_t q) { static const float C = 0x9.81b88p-4f; static const float P[] = {0x9.f2a54p-4f, 0x2.2b9cfp-4f, 0x7.32ff48p-12f}; static const float Q[] = {0x7.9a5a8p-4f, 0xe.7215p-8f}; float w = fabsf(q[3]); float n = P[2]; float d = Q[1]; n = fmaf(n,w,P[1]); n = fmaf(n,w,P[0]); n = fmaf(n,w,C); n = copysignf(n,q[3]); d = fmaf(d,w,Q[0]); d = fmaf(d,w,1.f); d = fmaf(d,w,C); return (n/d)*quat_bivector(q); }

연결된 파일에 적힌 대로, 위 유리 근사는 rminimax로 다음 명령을 사용해 만들었다.

ratapprox --function="2*(acos(x)/(sqrt(1-x^2)))/pi" --dom=[0,0.99999999] --denF=[SG] --numF=[SG] --num=[1,x,x^2,x^3] --den=[1,x,x^2,x^3] --weight=1 --output=fem_33.sollya

지수 공간에서 선형 공간으로

특수화된 역변환은 $\mathbf{v} = \theta~\mathbf{u}$일 때 다음과 같다.

[\begin{align*}

Q & = \tfunc{exp}{\frac{\pi}{2}~\mathbf{v}} \

& = \tfunc{cos}{\frac{\pi}{2}~\theta} + \tfunc{sin}{\frac{\pi}{2}~\theta}~\mathbf{u} \

& = \tfunc{cos}{\frac{\pi}{2}~\theta} + \frac{\tfunc{sin}{\frac{\pi}{2}~\theta}}{\sqrt{\mathbf{v} \cdot \mathbf{v}}}~\mathbf{v}

\end{align*}]

표준 함수를 이용한 기본 구현은 다음과 같다.

// exp(π/2 V) special cased for 'V' in unit ball // given V = ΘU (in unit ball, U = unit bivector) // returns cos(π/2 Θ) + sin(π/2 Θ) U // distance error: ~2.63713986e-07 static inline quatf_t quatf_iem_naive(vec3f_t v) { float d = vec3_norm_fma(v); // v∙v float x = sqrtf(d); // |v| float a = ((float)(0.5*M_PI))*x; // this casted double happens to be correctly rounded in singles float w = cosf(a); // or sincospi if you have it float s = sinf(a)/x; s = (x != 0.f) ? s : 1.f; // correct any limit case return quatf_bs(s*v,w); }

직접 구현하기 위해 $\cos$가 짝함수라는 점을 이용해 다음과 같이 근사한다.

[\tfunc{cos}{\frac{\pi}{2}~\theta} = 1 + \theta^2~\func{P}{\theta^2}]

여기서 $P$는 다항식이다. $\theta$를 계산하는 첫 단계는 다음과 같다.

[d = \mathbf{v} \cdot \mathbf{v} = \theta^2]

이를 코사인 식에 대입하고 $\theta$가 양수임을 이용해 $p = \func{P}{d}$로 놓으면 다음을 얻는다.

[\begin{align*}

\tfunc{cos}{\frac{\pi}{2}~\theta} & = 1 + dp \

\tfunc{sin}{\frac{\pi}{2}~\theta} & = \sqrt{1-\left(1 + dp\right)^2} \

& = \sqrt{-2dp-d^2~p^2}

\end{align*}]

따라서 다음 식이 성립한다.

[\begin{align*}

\frac{\tfunc{sin}{\frac{\pi}{2}~\theta}}{\sqrt{d}} & = \frac{\sqrt{-2dp-d^2~p^2}}{\sqrt{d}} \

& = \sqrt{\frac{-2dp-d^2~p^2}{d}} \

& = \sqrt{-2p-d~p^2}

\end{align*}]

정방향 함수와 마찬가지로 파일에 여러 구현을 담았다. 현재 세 가지가 있으며, 여기서는 앞서 본 구현보다 오차 한계가 작은 최소한의 예를 보인다.

// ~2.30486819e-07 static inline quatf_t quatf_iem(vec3f_t v) { // cos(π/2 Θ) ≈ 1 + Θ²P(Θ²) where Θ² = d = v∙v // coefficients of polynomial P static const float C[] = {0x1.c2a574p-11f, -0x1.5502d4p-6f, 0x1.03bd9ep-2f, -0x1.3bd3bp0f}; // computation of 's' // sin(π/2 Θ)/Θ = sqrt(1-w²)/Θ : from cosine // = sqrt((1-w²)/d) : pull Θ inside sqrt // = sqrt((1-(1+dp)²)/d) : w = 1+dp // = sqrt(-dp²-2p) : reduce float d = vec3_norm_fma(v); // v∙v float p = f32_horner_3(d,C); // p = P(d) float s = sqrtf(-fmaf(d,p*p,p+p)); // sin(π/2 Θ)/Θ * float w = fmaf(d,p,1.f); // cos(π/2 Θ) * // { sin(π/2 Θ)/Θ v, cos(π/2 Θ) } return quatf_bs(s*v,w); }

별표가 붙은 두 줄은 오차를 줄이도록 FMA 연산을 구성했다. 그래서 $pd$를 두 번 계산한다. FMA 없이 이 값들을 계산하는 가장 좋은 방법은 따로 검토하지 않았다.

결과

균일한 무작위 단위 쿼터니언을 만들고, 지수 공간으로 변환한 뒤 다시 선형 공간으로 되돌려 원래 쿼터니언과의 3차원 회전각을 측정했다. 회전각을 라디안으로 계산하는 코드는 다음과 같다.

// computes the 3D rotation (instead of intrinsic) angle between // A & B measured in radians. static inline double quatf_relative_angle_hq(quatf_t A, quatf_t B) { quatd_t a = quatf_promote(A); // promote to doubles quatd_t b = quatf_promote(B); // ditto quatd_t r = quatd_mul(a,quat_conj(b)); // compute relative: AB^* // use atan2 for higher accuracy return 2.0*atan2(sqrt(vec3_dot_fma(r,r)),r[3]); }

다음 정방향·역방향 함수 조합을 비교했다.

• 기준 구현: 연결된 파일에 있는 표준 함수를 사용하며, 배정밀도로 변환해 계산한다.
• 교란: 변환 없이 네 성분 각각의 최하위 비트만 뒤집는다.
• 표준 함수 구현.
• 이 글에서 선택한 직접 구현.

기준 함수 조합은 각 성분의 비트를 교란하는 것보다 오차가 작다. 예시로 든 직접 구현은 표준 함수 구현보다 약간 나쁘지만, 절대 차이는 $9 \times 10^{-9}$보다 작다. 이 결과를 성공이라고 볼 만하다.

쿼터니언 관련 글

애니메이션, 회전