core 設計

このページでは、Luna-Flow/quaternion が実装する数学と、API が現在の形になっている理由を説明します。以下の式はすべて src/quaternion.mbt のコードが実際に評価する式であり、実装が数学から外れている箇所はそのことを明記します。

設計目標

本パッケージは MoonBit に四元数型 Quaternion[T] を 1 つ提供し、2 種類の利用者に応えます。

• 代数的なコード：$\mathbb H$ を厳密または近似のスカラー型上の環として扱い、luna-generic の trait を通して利用するもの；
• 幾何的なコード：単位四元数を 3 次元回転として扱い、軸と角またはオイラー角から作成し、合成・適用・補間・分解するもの。

どちらも同じジェネリック型で扱えます。各演算は必要な T の trait だけを要求するため、環演算は Int 上で厳密に行え、平方根や三角関数が必要な演算は Double を経由します。

数学的背景

実代数としての四元数

四元数 $\mathbb H$ は、$1, i, j, k$ を基底とする 4 次元実ベクトル空間で、

$$q = w + x\,i + y\,j + z\,k, \qquad w, x, y, z \in \mathbb R,$$

その双線形な乗法は Hamilton の関係式で定まります¹1 W. R. Hamilton, “On quaternions; or on a new system of imaginaries in algebra”, 1843 年。これらの関係式は、彼がダブリンのブルーム橋（Brougham Bridge）に刻んだものです。

$$i^2 = j^2 = k^2 = ijk = -1 .$$

$q = (w, \mathbf u)$ と書き、$w$ をスカラー部、$\mathbf u = (x, y, z)$ をベクトル部とし、$\mathbb R^3$ を純四元数 $\{(0, \mathbf u)\}$ と同一視します。Quaternion[T] はまさにこの組を保持します：$w$ にはフィールド r、$\mathbf u$ には三つ組 vec を使います。乗法は $\mathbb R$-双線形なので、実数 $(s, \mathbf 0)$ はすべての元と可換です。

基底同士の積

積の性質はすべてこの 4 つの関係式から導かれます。$ijk = -1$ の右から $k$ を掛け、$k^2 = -1$ を使います。

$$\begin{aligned} ijk\,k &= -k \\ -ij &= -k \\ ij &= k . \end{aligned}$$

$ijk = -1$ の左から $i$ を掛け、$i^2 = -1$ を使うと $-jk = -i$、よって $jk = i$ です。逆順の積はこの 2 つから得られます。

$$\begin{aligned} ji &= j\,(jk) = j^2 k = -k, \\ kj &= (ij)\,j = i\,j^2 = -i, \\ ik &= i\,(ij) = i^2 j = -j, \\ ki &= (ij)\,i = i\,(ji) = -ik = j . \end{aligned}$$

したがって基底同士の積は、標準基底の外積に対角成分 $-1$ を加えたものになります。

$\cdot$	$i$	$j$	$k$
$i$	$-1$	$k$	$-j$
$j$	$-k$	$-1$	$i$
$k$	$j$	$-i$	$-1$

8 つの元 $\pm 1, \pm i, \pm j, \pm k$ は四元数群 $Q_8$ をなし、$2\times 2$ 複素行列で実現できます。²2 例えば $1 \mapsto I$、$i \mapsto \begin{pmatrix} i & 0 \\ 0 & -i\end{pmatrix}$、$j \mapsto \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix}$、$k \mapsto \begin{pmatrix} 0 & i \\ i & 0\end{pmatrix}$。これらの行列は Hamilton の関係式を満たし、行列の乗法は結合的です。 $Q_8$ の乗法は結合的であり、双線形性から $\mathbb H$ の積も結合的です。

Hamilton 積

表を使って $p = w_1 + x_1 i + y_1 j + z_1 k$ と $q = w_2 + x_2 i + y_2 j + z_2 k$ の積を項ごとに展開すると、次が得られます。

$$\begin{aligned} pq ={} & (w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2) \\ &+ (w_1 x_2 + x_1 w_2 + y_1 z_2 - z_1 y_2)\, i \\ &+ (w_1 y_2 - x_1 z_2 + y_1 w_2 + z_1 x_2)\, j \\ &+ (w_1 z_2 + x_1 y_2 - y_1 x_2 + z_1 w_2)\, k . \end{aligned}$$

同じ積はスカラー–ベクトル形式の方が見通しがよくなります。2 つの純四元数 $\mathbf a$ と $\mathbf b$ について、基底の二乗は $-a_1 b_1 - a_2 b_2 - a_3 b_3$ を与え、混合項は対になります。例えば $a_1 b_2\, ij + a_2 b_1\, ji = (a_1 b_2 - a_2 b_1)k$ です。したがって

$$\mathbf a\,\mathbf b = -\,\mathbf a\cdot\mathbf b + \mathbf a\times\mathbf b .$$

双線形性とスカラーが可換であることから、

$$\begin{aligned} (s_1, \mathbf v_1)(s_2, \mathbf v_2) &= s_1 s_2 + s_1\mathbf v_2 + s_2\mathbf v_1 + \mathbf v_1\mathbf v_2 \\ &= \big(s_1 s_2 - \mathbf v_1\cdot\mathbf v_2,\ \ s_1\mathbf v_2 + s_2\mathbf v_1 + \mathbf v_1\times\mathbf v_2\big). \end{aligned}$$

これはまさに Mul の実装です。スカラー部は self.r * other.r - dot(self.vec, other.vec)、ベクトル部は cross(self.vec, other.vec) に 2 つのスケーリングしたベクトルを加えたものです。

因子を入れ替えると外積の符号だけが反転するので、

$$pq - qp = 2\,\mathbf v_1 \times \mathbf v_2 ,$$

2 つの四元数が可換であるのは、ベクトル部が平行なときに限ります。特に $\mathbb H$ は可換ではありません：$ij = k$ ですが $ji = -k$ です。

この導出では成分同士が可換であること（$s_1\mathbf v_2 = \mathbf v_2 s_1$ など）を使いました。そのためジェネリックな実装は T の乗法が可換であることを前提とし、luna-generic が提供するすべての数値型はこれを満たします。

共役とノルム

共役は $\bar q = (w, -\mathbf u)$ です。平行なベクトルの外積は零なので、四元数にその共役を掛けると実数が残ります。

$$q\bar q = \big(w^2 + \mathbf u\cdot\mathbf u,\ -w\mathbf u + w\mathbf u - \mathbf u\times\mathbf u\big) = \big(w^2 + x^2 + y^2 + z^2,\ \mathbf 0\big) = |q|^2 ,$$

同様に $\bar q q = |q|^2$ です。Quaternion::square_len はこの $|q|^2$ を計算し、Quaternion::dot はそれをノルムとする $\mathbb R^4$ の内積です。

共役は積の順序を反転します。$p = (s_1, \mathbf v_1)$、$q = (s_2, \mathbf v_2)$ とすると、

$$\begin{aligned} \overline{pq} &= \big(s_1 s_2 - \mathbf v_1\cdot\mathbf v_2,\ -s_1\mathbf v_2 - s_2\mathbf v_1 - \mathbf v_1\times\mathbf v_2\big), \\ \bar q\,\bar p &= \big(s_2 s_1 - \mathbf v_2\cdot\mathbf v_1,\ -s_2\mathbf v_1 - s_1\mathbf v_2 + \mathbf v_2\times\mathbf v_1\big), \end{aligned}$$

$\mathbf v_2\times\mathbf v_1 = -\mathbf v_1\times\mathbf v_2$ なので両者は一致します。よって $\overline{pq} = \bar q\,\bar p$ です。

ノルムは乗法的です。結合性、$\overline{pq} = \bar q\bar p$、および実数 $|q|^2$ が $p$ と可換であることを使うと、

$$\begin{aligned} |pq|^2 &= (pq)\,\overline{(pq)} = p\,(q\bar q)\,\bar p \\ &= p\,|q|^2\,\bar p = |q|^2\,p\bar p = |p|^2\,|q|^2 . \end{aligned}$$

整数上ではこれはオイラーの四平方恒等式であり、チュートリアルで Quaternion[Int] を使って確かめています。

逆元と 2 つの除算

$q$ が零でなければ $|q|^2 > 0$ であり、$q\bar q = \bar q q = |q|^2$ より

$$q^{-1} = \frac{\bar q}{|q|^2}$$

は両側逆元です。これが Inverse の実装で、共役を one() / square_len() でスケーリングします。零でない四元数はすべて可逆なので、$\mathbb H$ は斜体（除環）です。可換でないため、体ではありません。

可換性がないと、「$q$ を $r$ で割る」には未知数がどちら側にあるかに応じて 2 つの意味があります。

$$\begin{aligned} x\,r = q &\iff x = q\,r^{-1} && \text{(right division, } q / r\text{)}, \\ r\,x = q &\iff x = r^{-1} q && \text{(left division, } q\text{.left\_div}(r)\text{)}. \end{aligned}$$

実装は $r^{-1}$ を作らず、最後に 1 回だけ割ります。$q = (s_q, \mathbf v_q)$、$r = (s_r, \mathbf v_r)$、$\bar r = (s_r, -\mathbf v_r)$ とすると、スカラー–ベクトル積から

$$\begin{aligned} q\,r^{-1} &= \frac{q\,\bar r}{|r|^2} = \frac{\big(s_q s_r + \mathbf v_q\cdot\mathbf v_r,\ \ s_r\mathbf v_q - s_q\mathbf v_r - \mathbf v_q\times\mathbf v_r\big)}{|r|^2}, \\ r^{-1} q &= \frac{\bar r\,q}{|r|^2} = \frac{\big(s_q s_r + \mathbf v_q\cdot\mathbf v_r,\ \ s_r\mathbf v_q - s_q\mathbf v_r - \mathbf v_r\times\mathbf v_q\big)}{|r|^2}, \end{aligned}$$

これらはそれぞれ Div の実装（c = cross(q.vec, r.vec)）と Quaternion::left_div（c = cross(r.vec, self.vec)）です。2 つの商はスカラー部が等しく、差は次のとおりです。

$$q\,r^{-1} - r^{-1} q = -\frac{2\,\mathbf v_q\times\mathbf v_r}{|r|^2}.$$

例。 $q = 1 + 2i + 3j + 4k$、$r = 5 + 6i + 7j + 8k$ とします。すると $|r|^2 = 174$、$s_q s_r + \mathbf v_q\cdot\mathbf v_r = 5 + 12 + 21 + 32 = 70$、$s_r\mathbf v_q - s_q\mathbf v_r = (4, 8, 12)$、$\mathbf v_q\times\mathbf v_r = (-4, 8, -4)$ なので、

$$\begin{aligned} q / r &= \tfrac{1}{174}\big(70 + 8i + 0j + 16k\big) \approx 0.4023 + 0.0460\,i + 0.0920\,k, \\ q.\mathrm{left\_div}(r) &= \tfrac{1}{174}\big(70 + 0i + 16j + 8k\big) \approx 0.4023 + 0.0920\,j + 0.0460\,k . \end{aligned}$$

API ページで、これらの値と 2 つの定義式を確かめています。

単位四元数と回転

単位四元数（$|q| = 1$）は次のように書けます。

$$q = \Big(\cos\tfrac{\theta}{2},\ \sin\tfrac{\theta}{2}\,\mathbf n\Big), \qquad |\mathbf n| = 1 ,$$

これはまさに from_axis_angle(axis, angle) が返すもので、$\mathbf n$ は正規化した軸です。単位 $q$ では $q^{-1} = \bar q$ です。純四元数上の写像 $\mathbf v \mapsto q\,\mathbf v\,\bar q$ を考えます。$q = (w, \mathbf u)$ と書きます。まず、

$$q\,\mathbf v = \big(-\mathbf u\cdot\mathbf v,\ w\mathbf v + \mathbf u\times\mathbf v\big).$$

右から $\bar q = (w, -\mathbf u)$ を掛けると、スカラー部は

$$-w\,\mathbf u\cdot\mathbf v + (w\mathbf v + \mathbf u\times\mathbf v)\cdot\mathbf u = -w\,\mathbf u\cdot\mathbf v + w\,\mathbf v\cdot\mathbf u + 0 = 0 ,$$

となるので像は再び純四元数であり、ベクトル部は

$$\begin{aligned} \mathbf v' &= (\mathbf u\cdot\mathbf v)\,\mathbf u + w\,(w\mathbf v + \mathbf u\times\mathbf v) - (w\mathbf v + \mathbf u\times\mathbf v)\times\mathbf u \\ &= (\mathbf u\cdot\mathbf v)\,\mathbf u + w^2\mathbf v + 2w\,\mathbf u\times\mathbf v + \mathbf u\times(\mathbf u\times\mathbf v) \\ &= (w^2 - |\mathbf u|^2)\,\mathbf v + 2(\mathbf u\cdot\mathbf v)\,\mathbf u + 2w\,\mathbf u\times\mathbf v , \end{aligned}$$

ここで $\mathbf a\times\mathbf u = -\mathbf u\times\mathbf a$ と $\mathbf u\times(\mathbf u\times\mathbf v) = (\mathbf u\cdot\mathbf v)\mathbf u - |\mathbf u|^2\mathbf v$ を使いました。$w = \cos\frac\theta2$、$\mathbf u = \sin\frac\theta2\,\mathbf n$ を代入し、半角の公式 $\cos^2\frac\theta2 - \sin^2\frac\theta2 = \cos\theta$、$2\sin^2\frac\theta2 = 1 - \cos\theta$、$2\sin\frac\theta2\cos\frac\theta2 = \sin\theta$ を使うと、

$$\mathbf v' = \cos\theta\,\mathbf v + (1 - \cos\theta)(\mathbf n\cdot\mathbf v)\,\mathbf n + \sin\theta\,(\mathbf n\times\mathbf v).$$

これは Rodrigues の回転公式です：$\mathbf v \mapsto q\mathbf v q^{-1}$ は $\mathbf n$ の周りの角 $\theta$ の回転で、右手の法則で反時計回りです。$|q\mathbf v\bar q| = |q||\mathbf v||\bar q| = |\mathbf v|$ なので、この写像は回転にふさわしく等長写像です。

ここから API を形作る 3 つの帰結が得られます。

• 合成は乗法。 $p\,(q\mathbf v q^{-1})\,p^{-1} = (pq)\,\mathbf v\,(pq)^{-1}$ なので、「先に $q$、次に $p$」の回転は p * q です。
• 二重被覆。 $(-q)\,\mathbf v\,(-q)^{-1} = q\mathbf v q^{-1}$ なので、$q$ と $-q$ は同じ回転です。実際、角 $\theta + 2\pi$ は $-q$ を与えます。各回転にはちょうど 2 つの単位四元数が対応します。== で回転を比較できず、slerp が一方の入力の符号を反転するのはこのためです。
• 安価な評価。 Quaternion::rotate は $\mathbf t = 2\,\mathbf u\times\mathbf v$ と $\mathbf v' = \mathbf v + w\,\mathbf t + \mathbf u\times\mathbf t$ を評価します。展開すると $\mathbf u\times\mathbf t = 2(\mathbf u\cdot\mathbf v)\mathbf u - 2|\mathbf u|^2\mathbf v$ なので、 $$\mathbf v + w\,\mathbf t + \mathbf u\times\mathbf t = (1 - 2|\mathbf u|^2)\,\mathbf v + 2(\mathbf u\cdot\mathbf v)\,\mathbf u + 2w\,\mathbf u\times\mathbf v ,$$ これが上の式と一致するのは $1 - 2|\mathbf u|^2 = w^2 - |\mathbf u|^2$、すなわち $w^2 + |\mathbf u|^2 = 1$ のときに限ります。rotate が単位四元数を要求するのはこのためです。除算は不要ですが、ノルムが 1 であることに依存しています。

同じ式を基底ベクトルに適用すると、単位四元数の回転行列が得られます。

$$R(q) = \begin{pmatrix} 1 - 2(y^2 + z^2) & 2(xy - wz) & 2(xz + wy) \\ 2(xy + wz) & 1 - 2(x^2 + z^2) & 2(yz - wx) \\ 2(xz - wy) & 2(yz + wx) & 1 - 2(x^2 + y^2) \end{pmatrix}.$$

本パッケージはこの行列を返しませんが、後述のオイラー角の抽出はその成分を読み取ります。

オイラー角

座標軸の周りに $\alpha$ だけ回転する四元数を $q_X(\alpha)$、$q_Y(\alpha)$、$q_Z(\alpha)$、その行列を $R_X, R_Y, R_Z$ と書きます。軸 $A, B, C$ の周りの角 $a, b, c$ の 3 回の回転の列は、2 通りに読めます。

• 外因的（external=true）：固定軸の周りに、$A$、$B$、$C$ の順に回転します。後の回転は左から掛かります：$q = q_C(c)\,q_B(b)\,q_A(a)$。
• 内因的（external=false）：物体とともに動く軸の周りに、$A$、回転後の $B$、2 回回転した後の $C$ の順に回転します：$q = q_A(a)\,q_B(b)\,q_C(c)$。

同じ積を 2 通りに読むと、角 $(a, b, c)$ の内因的 $ABC$ は角 $(c, b, a)$ の外因的 $CBA$ であることがわかります。実装はまさにこれを使っており、各 to_euler_internal_* 関数は逆順の外因的関数を呼び、三つ組を反転します。

抽出：外因的 XYZ。 $R = R_Z(c)\,R_Y(b)\,R_X(a)$ について、基本行列を掛け合わせると

$$R = \begin{pmatrix} \cos b\cos c & \cdots & \cdots \\ \cos b\sin c & \cdots & \cdots \\ -\sin b & \cos b\sin a & \cos b\cos a \end{pmatrix}.$$

$R(q)$ と成分ごとに比較すると、

$$\begin{aligned} \sin b &= -R_{31} = 2(wy - xz), \\ a &= \operatorname{atan2}(R_{32}, R_{33}) = \operatorname{atan2}\big(2(wx + yz),\ 1 - 2(x^2 + y^2)\big), \\ c &= \operatorname{atan2}(R_{21}, R_{11}) = \operatorname{atan2}\big(2(wz + xy),\ 1 - 2(y^2 + z^2)\big), \end{aligned}$$

これは $\cos b > 0$ の間成り立ちます（atan2 の両引数を $\cos b$ で割っても角は変わりません）。これらが to_euler_external_XYZ の式で、$b = \arcsin(\cdot) \in [-\pi/2, \pi/2]$ です。to_euler_external_YZX と to_euler_external_ZYX は軸を置換した同じ導出です。巡回的でない置換では、逆正弦の中身と atan2 の分子の符号が反転します。モジュールは抽出前に q を正規化するので、$R(q)$ の対角成分で使った恒等式 $w^2 + x^2 + y^2 + z^2 = 1$ は丸め誤差の範囲で成り立ちます。

ジンバルロック。 $\cos b = 0$ のとき、成分 $R_{32}, R_{33}, R_{21}, R_{11}$ はすべて零になり、1 回目と 3 回目の回転は同じ物理的な軸の周りに作用します。XYZ の場合、$b = \pi/2$ では

$$R_Z(c)\,R_Y(\tfrac\pi2)\,R_X(a) = R_Y(\tfrac\pi2)\,R_X(a - c),$$

となり、決まるのは $a - c$ だけです（$b = -\pi/2$ では $a + c$）。慣例的な対処は $c = 0$ とし、残った成分から $a$ を求めることで、ここでは $a = \operatorname{atan2}(R_{12}, R_{22})$ です。実装は $|\sin b| \ge 0.9998$（約 $88.85^\circ$）で特別な分岐に切り替わり、警告を出力して $b = \pm\pi/2$、$c = 0$ としますが、最初の角は XYZ では $\operatorname{atan2}(R_{13}, R_{22})$ から、その他の順序では通常の分岐の式（引数がほぼ零）から求めます。どちらも $a \mp c$ を分離しないため、この分岐は現在、入力を再現しない角を返します。既知の逸脱を参照してください。

構築。 既定の "XYZ" での from_euler(roll, pitch, yaw) は、$q_X(\text{roll})\,q_Y(\text{pitch})\,q_Z(\text{yaw})$ を閉じた形で展開します。$c_r = \cos\frac{\text{roll}}{2}$、$s_r = \sin\frac{\text{roll}}{2}$ などとすると、Hamilton 積を 2 回適用して

$$q_X q_Y q_Z = \big(c_r c_p c_y - s_r s_p s_y,\ \ s_r c_p c_y + c_r s_p s_y,\ \ c_r s_p c_y - s_r c_p s_y,\ \ c_r c_p s_y + s_r s_p c_y\big),$$

これがソースの式です。to_euler_internal_XYZ はその逆です。

球面線形補間

単位四元数は 3 次元球面 $S^3 \subset \mathbb R^4$ をなし、2 つの姿勢の間の最短経路は大円の弧です。$q_1, q_2$ を単位四元数、$\cos\Omega = q_1\cdot q_2$、$0 < \Omega < \pi$ とします。単位ベクトル

$$q_\perp = \frac{q_2 - \cos\Omega\, q_1}{\sin\Omega}$$

は $q_1$ と直交し、弧は $\gamma(\varphi) = \cos\varphi\, q_1 + \sin\varphi\, q_\perp$ です。$\varphi = t\Omega$ では、

$$\begin{aligned} \gamma(t\Omega) &= \frac{\sin\Omega\cos t\Omega - \cos\Omega\sin t\Omega}{\sin\Omega}\,q_1 + \frac{\sin t\Omega}{\sin\Omega}\,q_2 \\ &= \frac{\sin\big((1-t)\Omega\big)}{\sin\Omega}\,q_1 + \frac{\sin(t\Omega)}{\sin\Omega}\,q_2 , \end{aligned}$$

これが slerp の評価する式です。一定の角速度で動き、回転角は $t$ に比例して増えます。二重被覆のため、slerp は $q_1\cdot q_2 < 0$ のときまず $q_2$ を $-q_2$ に置き換えて $\Omega \le \pi/2$ とし、回転が短い方を回るようにします。丸めによって 2 つの単位四元数の内積がわずかに 1 を超えることがあるため、acos の前に $\cos\Omega$ を $[-1, 1]$ に制限します。

べき乗

pow_by_int は $q^{2m} = (q^m)^2$ と $q^{2m+1} = q\,(q^m)^2$ を使い、結合性だけを必要とします。すべての因子は同じ $q$ のべきなので互いに可換です。負の指数には $q^{-n} = (q^{-1})^n$ を使います。

pow_by_T は極形式を使います。零でない $q$ はすべて $q = |q|\,(\cos\varphi + \hat{\mathbf n}\sin\varphi)$（$\varphi \in [0, \pi]$、$\hat{\mathbf n}$ は単位ベクトル）と書け、$1$ と $\hat{\mathbf n}$ が張る空間は $\mathbb C$ の複製です（$\hat{\mathbf n}^2 = -1$ のため）。するとド・モアブルの公式により次が定義されます。

$$q^t = |q|^t\,\big(\cos t\varphi + \hat{\mathbf n}\sin t\varphi\big).$$

実装は $\varphi = \arcsin\lvert\mathbf u\rvert/|q|$ を計算しますが、これが正しいのは $\varphi \le \pi/2$、すなわちスカラー部が非負の場合だけです。既知の逸脱を参照してください。

設計上の決定

ジェネリックな成分と関数ごとの制約

問題。 四元数は厳密な整数上（数論、恒等式の検証）でも Double 上（幾何）でも有用であり、luna-generic のエコシステムは両方に合う 1 つの型を必要としています。

選択肢。 Double 専用の型、単一の「実数」trait でパラメータ化した型、あるいは T に制約を課さず各関数が使うものだけを要求する型。

選択。 最後のものです。Quaternion[T] には制約がなく、+ は Add、* は Mul + Sub + Add、square_len は Mul + Add を要求し、平方根や三角関数を使う関数だけが DoubleConvert を要求します。これは、必要なものを表す最小の trait の組み合わせに依存するという Luna-Flow の原則に従ったもので、Quaternion[Int] で代数的恒等式を厳密に確かめることもできます。

スカラー–ベクトルによる格納

この型は 4 つの独立したフィールドではなく、r : T と vec : (T, T, T) を保持します。上の積・共役・逆元・回転はスカラー–ベクトル形式で導出されているので、コードは導出そのもののように読めます（vec に対する dot、cross、scale）。型は抽象型なので表現を自由に変えられますが、その代償として成分のアクセサがまだありません。

Field ではなく Ring である理由

luna-generic の Field は Ring + Inverse + Div であり、それに対して書かれたジェネリックなコードは、$a b = b a$ や $a / b \cdot c = a c / b$ のような可換体の法則に頼ってよいことになっています。四元数はすべての環の公理を満たし逆元も持ちますが、可換ではないため、Field を名乗るとそのようなコードが黙って誤った答えを出しかねません。そこで本パッケージは Zero、One、AddMonoid、MulMonoid、Semiring、Ring（可換な T 上の $\mathbb H$ ではすべて成り立つ）に加えて演算 trait の Inverse と Conjugate を実装し、そこで止めています。除算は Div と left_div で引き続き使えますが、ジェネリックなコードは明示的に要求する必要があります。

/ は右除算

問題。 乗法が可換でない場合、q / r はどちらかの側を選ぶ必要があります。0.2.0 より前は $r^{-1} q$ を計算していました。

選択。 0.2.0 から q / r は $q\,r^{-1}$、すなわち $x r = q$ の解です。これは斜体での通常の約束であり、luna-generic のスカラー型と同じように Div と Inverse の整合性を保ちます：a / b == a * b.inv()。回転としても自然に読めます：to / from は from の後に適用すると to になる回転です。もう一方の商は、逆元がどちら側にあるかを名前で示す left_div として引き続き提供されます。この変更は破壊的変更であり、変更履歴に記録されています。

Double への橋渡し

平方根と三角関数は Double にしかなく、本パッケージは解析的な trait を持たない luna-generic にしか依存していません。DoubleConvert は成分を Double に変換して戻す 2 メソッドの trait で、これにより magnitude、normalize、slerp、pow_by_T、オイラー角変換は、Double で計算しながらシグネチャ上はジェネリックなままでいられます。Double（恒等写像）と Int（from_double は切り捨て）に実装されています。この trait は pub(open) ではなく pub として宣言されているため、他のパッケージが自分の型に実装することはできません。

厳密な構造的等価性

Eq と Hash は導出されたもので、成分を厳密に比較します。「回転として」の等価性（$q \sim -q$）や「許容誤差内」の等価性は用途によって異なるため、呼び出し側に委ねます（チュートリアルで 1 つの方法を示しています）。厳密な等価性により Eq と Hash の整合性も保たれます。

明示的なメソッドの昇格

MoonBit 0.10 では trait の実装が暗黙にメソッドになることはなくなりました。src/extends.mbt は、四元数の自然な演算である演算子、equal、hash、to_string、conjugate、zero、one、inv を昇格します。メソッド形式の not_equal、hash_combine、output、default、to_repr は互換性のために残され、非推奨として隠されています。

正しさと不変条件

代数的な法則

厳密な算術を持つ可換環 T（Int、BigInt）上で、実装は次を満たします。

• 環の公理：$(\mathbb H, +, 0)$ はアーベル群、$(\mathbb H, \cdot, 1)$ はモノイドであり、$\cdot$ は $+$ に対して両側分配的；
• $\overline{pq} = \bar q\,\bar p$、$\bar{\bar q} = q$、$q\bar q = \bar q q = |q|^2$、$|pq|^2 = |p|^2|q|^2$；
• 厳密な除算（有理数のように T が体の場合）：$q\,q^{-1} = q^{-1}q = 1$、$(q / r)\,r = q$、$r\,(q.\mathrm{left\_div}(r)) = q$。

パッケージのテストは、整数のサンプルで環の公理・非可換性・2 つの除算の恒等式を確かめています。

浮動小数点誤差

Double では、これらの法則は丸め誤差の範囲で成り立ちます。$pq$ の各成分は、$p$ の成分と $q$ の成分の 4 つの積の和であり、積で 4 回、加算で 3 回の丸めを伴って評価されます。このような内積に対する標準的な誤差限界³3 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 第 2 版, SIAM 2002, §3.1：任意の和の順序で $\lvert\mathrm{fl}(x^{\mathsf T}y) - x^{\mathsf T}y\rvert \le \gamma_n \lvert x\rvert^{\mathsf T}\lvert y\rvert$。 は、成分ごとに次を与えます。

$$\lvert \mathrm{fl}(pq)_m - (pq)_m \rvert \le \gamma_4 \sum_{l=1}^{4} \lvert p_l\rvert\,\lvert q_{\sigma_m(l)}\rvert \le \gamma_4\, |p|\,|q|, \qquad \gamma_n = \frac{n u}{1 - n u},$$

ここで $\sigma_m$ は成分 $m$ に現れる $q$ の成分の置換であり、2 つ目の不等式は Cauchy–Schwarz の不等式です。4 成分について和をとると $\lVert\mathrm{fl}(pq) - pq\rVert \le 2\gamma_4 |p||q|$ となり、ノルムの乗法性から

$$(1 - 2\gamma_4)\,|p|\,|q| \le |\mathrm{fl}(pq)| \le (1 + 2\gamma_4)\,|p|\,|q| .$$

したがって $n$ 個の単位四元数の積のノルムは $1 + \delta$ で、一次近似で $|\delta| \lesssim 8nu$ です（$u = 2^{-53} \approx 1.1 \times 10^{-16}$）。実際には誤差が部分的に打ち消し合い、ずれはこれより小さくなります。rotate は $|q| = 1$ を前提とするため、ずれた四元数はベクトルを約 $1 + 2\delta$ 倍に拡大縮小します。長い積の連鎖では定期的に normalize を呼ぶべきです。その結果のノルムは $u$ の数倍以内で $1$ になります（例えば $1 + 2i + 3j + 4k$ を正規化すると、計算上のノルムは $0.9999999999999999$ になります）。

ゼロ除算と退化した入力

本パッケージは入力を検査しません。退化した場合は T の算術に従います。

入力	Double	Int
除数が零の inv、/、left_div	NaN 成分	実行時トラップ（整数のゼロ除算）
零の normalize	零をそのまま返す	零をそのまま返す
零の pow_by_T	零をそのまま返す	零をそのまま返す
軸が零の from_axis_angle	NaN のベクトル部	実行時トラップ
入力が零の slerp	零の入力が正規化されずにそのまま使われる	意味を持たない

Int では $|q|^2 = 1$ でない限り inv は零になり、DoubleConvert を経由する関数はすべて切り捨てを行います。

しきい値

slerp は $\cos\Omega > 0.9995$、すなわち $\Omega < 0.0316$ rad のとき、正規化した線形補間に切り替わります。そこでは $\sin\Omega < 0.032$ であり、これで割ると重みの丸め誤差が増幅されます。一方、しきい値において線形な経路と弧のずれは $S^3$ 上で $5.1 \times 10^{-7}$ rad 未満（回転角にして約 $10^{-6}$ rad）であり、しきい値より下ではさらに小さくなります。

オイラー角の変換は $|\sin b| \ge 0.9998$ をジンバルロックとみなします。そこで $b$ を $\pm\pi/2$ に丸めると、中央の角に最大 $\pi/2 - \arcsin 0.9998 \approx 0.020$ rad の誤差が生じます。

計算量

すべての演算は 4 成分に対する $O(1)$ です。Hamilton 積は乗算 16 回と加算 12 回、rotate は乗算 18 回と加算 12 回です。pow_by_int(n) は $O(\log |n|)$ 回の積と同程度の再帰の深さを使います。

既知の逸脱

以下は、現在の実装が上記の数学から外れている点です。隠さずにここに記載し、修正のために報告済みです。

• "ZYX" または "YZX" の from_euler は軸回転の積を評価しておらず、結果は単位四元数になりません。"ZXY" は角を位置で、"YXZ" は軸で割り当てます。
• to_euler_external_XZY（したがって to_euler_internal_YZX も）は、指定した系列ではなく内因的 X-Y-Z の抽出を評価します。
• オイラー角の抽出におけるジンバルロックの分岐は自由な角を分離せず（上記参照）、警告を println で出力します。
• pow_by_T は $\varphi$ を asin で計算するため $\pi/2$ を超えられません。$\operatorname{atan2}(\lvert\mathbf u\rvert, w)$ を使えば $[0, \pi]$ 全体を扱えます。

採用しなかった案

• Field の実装。 ジェネリックな Field のコードが可換性に依存しうるため採用しませんでした（上記参照）。
• / を左除算のままにする。 a / b == a * b.inv() にも、回転における to / from の通常の読み方にも反していました。
• Double 専用の型。 Int 上の厳密な算術と、他のスカラー型での luna-generic の環構造が失われます。
• 4 つの名前付きフィールドや回転行列での格納。 4 つのフィールドでは式が使うスカラー–ベクトル構造が隠れます。行列は 9 成分を持ち、安価な正規化の代わりに再直交化が必要で、補間も簡単ではありません。
• Eq での許容誤差による比較。 許容誤差に基づく == は推移的でなく、Hash と整合させることもできません。

境界

本パッケージは意図的に次のことを行いません。

• 回転行列、指数写像・対数写像、四元数の微積分（双対四元数、姿勢の微分）を提供すること；
• 入力を検査したり Result を返したりすること：退化した入力は T の算術に従い、rotate はレシーバが単位四元数であると信頼します；
• luna-generic の Field や Num を実装すること、あるいは四元数を順序付きの型にすること；
• arithmetic、linear-algebra、luna-complex に依存すること：唯一の Luna-Flow の依存先は luna-generic です；
• 回転の比較に使う許容誤差を決めること：それは呼び出し側に委ねます。

Footnotes

1. W. R. Hamilton, “On quaternions; or on a new system of imaginaries in algebra”, 1843 年。これらの関係式は、彼がダブリンのブルーム橋（Brougham Bridge）に刻んだものです。 ↩

2. 例えば $1 \mapsto I$、$i \mapsto \begin{pmatrix} i & 0 \\ 0 & -i\end{pmatrix}$、$j \mapsto \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix}$、$k \mapsto \begin{pmatrix} 0 & i \\ i & 0\end{pmatrix}$。これらの行列は Hamilton の関係式を満たし、行列の乗法は結合的です。 ↩

3. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 第 2 版, SIAM 2002, §3.1：任意の和の順序で $\lvert\mathrm{fl}(x^{\mathsf T}y) - x^{\mathsf T}y\rvert \le \gamma_n \lvert x\rvert^{\mathsf T}\lvert y\rvert$。 ↩