Rubyで四元数クラス（Quaternion）を実装するメモ

目的

• Rubyの数値クラスの仕組みを知る
  • 特に、異なるクラス間での演算はどう実現しているのか
  • 新しい数値クラスの導入にはどう備えているのか
• 四元数クラスを手に入れる
  • 要件として、組み込みの複素数クラス Complex と遜色なく扱えること
• 今後、4次元正多胞体で遊ぶ
  • 例：全頂点の座標を生成

完成品 → quaternion_c2（GitHub / RubyGems）

TL;DR：数値クラスの仕組み

自分のクラスで扱い方の分からないクラスのインスタンスが引数に来たら、そのインスタンスのメソッドを上手に呼び出して相手に処理してもらっている。これによって新しい数値クラスを追加してもうまく動くようになっている。

• 四則演算では、強制型変換 #coerce を使う
• Kernel.#Complex などコンストラクタで受け入れたいときは、四則演算に還元している
• 相手のクラス名を知らなくてもその性質は分かるように、 #real? など問い合わせメソッドがある

四元数の表現方法と計算式

インスタンス変数として何を記憶させると実装しやすいか考えるために、よくある表現方法（行列は除く）についての共役・絶対値・乗算の式をまとめておく。

実数4つ

基本となる方法。乗算がつらい…

\begin{align}
q &= w + xi + yj + zk \quad (w, x, y, z \in \mathbb{R}) \\
q^* &= w - xi - yj - zk \\
|q| &= \sqrt{w^2 + x^2 + y^2 + z^2} \\
q_1 q_2 &= (w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2) \\
& \;\;\;\; + (x_1 w_2 + w_1 x_2 - z_1 y_2 + y_1 z_2)i \\
& \;\;\;\; + (y_1 w_2 + z_1 x_2 + w_1 y_2 - x_1 z_2)j \\
& \;\;\;\; + (z_1 w_2 - y_1 x_2 + x_1 y_2 + w_1 z_2)k
\end{align}

複素数2つ

ケーリー＝ディクソンの構成法に基づく。複素数の計算に落とし込める。

\begin{align}
q &= a + bj \quad (a, b \in \mathbb{C}) \\
q^* &= a^* - bj \\
|q| &= \sqrt{|a|^2 + |b|^2} \\
q_1 q_2 &= (a_1 a_2 - b_2^* b_1) + (b_2 a_1 + b_1 a_2^*)j
\end{align}

この構成法を使えば 複素数→四元数→八元数→… とどんどん拡張することもできる。

スカラーと3次元ベクトル

3次元空間への応用が楽な形式。3次元ベクトルの計算に落とし込める。四元数が実部と虚部に分かれるため、極形式との変換でも必要になる。

\begin{align}
q &= s + \vec{v} \quad (s \in \mathbb{R}, \, \vec{v} \in \mathbb{R}^3) \\
q^* &= s - \vec{v} \\
|q| &= \sqrt{s^2 + |\vec{v}|^2} \\
q_1 q_2 &= (s_1 s_2 - \vec{v}_1 \cdot \vec{v}_2) + (s_1 \vec{v}_2 + s_2 \vec{v}_1 + \vec{v}_1 \times \vec{v}_2)
\end{align}

極形式（絶対値・角度・軸）

この表現のまま四則計算はできないが、角度や軸が出てくるため回転を考えるときに便利。また、累乗や指数・対数を計算するときに必要になる。

\begin{align}
q &= r \exp{(\vec{n} \theta)} \quad (r, \theta \in \mathbb{R}, \, \vec{n} \in \mathbb{R}^3, |\vec{n}| = 1) \\
  &= r \cos{\theta} + \vec{n} r \sin{\theta} \\
q^* &= r \exp{(-\vec{n} \theta)} \\
|q| &= |r| \\
q_1 q_2 & \neq r_1 r_2 \exp{(\vec{n}_1 \theta_1 + \vec{n}_2 \theta_2)} \\
q^\alpha &= r^\alpha \exp{(\vec{n} \alpha \theta)} \quad (\alpha \in \mathbb{R})
\end{align}

式の通り、一般に乗算を指数の加算に直すことはできない。乗算が可換、つまり $q_1 q_2 = q_2 q_1$ のときだけは等号が成立するらしい（参考）。

他の形式から変換するときは、一意に定まるように r≧0 や 0≦θ≦π と範囲を制限する。それでもqが0や負の実数のときは $\vec{n}\theta$ が無数に選べることに注意。

実装

今回の実装では複素数2つを記憶させる方法をとる。計算をなるべく Complex に還元することで実装の手間を省けると考えられる。

仕様についてはできる限り Complex を真似る。また、他の数値クラスとして BigDecimal や八元数を想定しておくと拡張性について考えやすい。

クラスの基盤

• 四元数クラス Quaternion は Numeric を継承
• インスタンス変数は複素数2つ
  • 二項演算を定義しやすいようgetterをprotectedで作成
  • 実数クラスの混在を避けるため、#initialize の中で必ず複素数に変換
• .new は外部から呼び出せないようprivateに設定

class Quaternion < Numeric
	attr_reader :a, :b
	protected   :a, :b

	def initialize(a, b)
		@a = a.to_c
		@b = b.to_c
	end

	private_class_method :new
end

定数：虚数単位

class Quaternion
	i = Complex::I # 一時変数
	I = new(i, 0)
	J = new(0, 1)
	K = new(0, i)
end

数の性質

• 体系
  • (#integer? ： 常に偽を返す)
  • #real? ： 常に偽を返す
  • #complex? [新規] ： 常に偽を返す
  • Numeric#complex? [新規] ： 常に真を返す
• 大きさ
  • #finite? ： 絶対値が有限かどうか返す
  • #infinite? ： 絶対値が無限なら 1、そうでなければ nil を返す
  • (#zero? ： self == 0)
  • (#nonzero? ： zero? ? nil : self)
  • #positive?, #negative?

例えば「あるオブジェクト num は実数と見なせるクラスのインスタンスか」を判定したいことがある。普通に考えると実数関連のクラスを列挙して [Integer, Rational, Float].any? { |c| num.kind_of?(c) } となるが、これでは新しい実数クラスが加わったときに対応できない。

そのためか数値クラスには問い合わせ用のメソッドがいくつかあり、実数の判定は num.kind_of?(Numeric) && num.real? とできる。新しい数値クラスを作るときはそれらのメソッドを適切に再定義する必要があり、四元数を作るのであれば #real? を偽にしておく。また、今回の実装では以降「複素数かどうか」が重要になるので、同様のメソッド #complex? を追加しておく。

実数かどうか判定

require 'bigdecimal'
num = BigDecimal('123.456')

[Integer, Rational, Float].any? { |c| num.kind_of?(c) } #=> false
num.kind_of?(Numeric) && num.real?                      #=> true

複素数は大小比較ができないが、大きさが有限か無限かについては判定が用意されている。あくまで絶対値に対する判定であるため、係数が有限かどうかは関係ない。四元数の場合も複素数への還元はせず、絶対値をうまく実装する必要がある。

有限かどうか判定

Complex(0         , 2 ** 1024 ).finite? #=> true
Complex(0.0       , 2 ** 1024 ).finite? #=> false
Complex(Float::MAX, Float::MAX).finite? #=> false

#zero? と #nonzero? は Numeric で定義されているので、#== さえ実装すれば正しく動作する。#positive? と #negative? は呼び出せないよう未定義にする。

単項演算

• (#+@ ： self)
• (#-@ ： 0 - self)
• #conj (#conjugate) ： 共役な数を返す
• #abs (#magnitude) ： 絶対値を返す
• #abs2 ： 絶対値の2乗を返す

#+@ と #-@ は Numeric で定義されている。処理を効率化するなら #-@ を再定義してもいい。

#abs と #abs2 は実装が独立している。#abs は Math.#hypot を使うことで、2乗すると浮動小数点数の範囲を超える場合でも計算できる。一方で #abs2 は単に実部と虚部の2乗を足すだけにして、数値の型変換を最小限に抑える。

代用不可な絶対値メソッド

c1 = Complex(Float::MAX, Float::MAX) / 2
c1.abs             #=> 1.2711610061536462e+308
Math.sqrt(c1.abs2) #=> Infinity

c2 = Complex(1, 2)
c2.abs2     #=> 5
c2.abs ** 2 #=> 5.000000000000001

四則演算

• #+, #-, #*, #/ ： 加減乗除¹
• #quo, #fdiv, #div ： 除算の種類
• #coerce ： 強制型変換

クラスが異なるもの同士で計算することはよくあるので、その状況に対応しなければならない。とはいえ決まったやり方に従えば割と簡単にできる。

まずは演算子の左側が自クラスの場合。メソッドの中で引数（演算子の右側）のクラスを見て処理を分ける。扱い方の分からないクラスだったときに、エラーは返さず相手の #coerce を呼ぶのがポイント。

乗算

class Quaternion
	def *(other)
		if other.kind_of?(Quaternion)
			# 四元数同士。前出の数式通りに計算すればいい。
			_a = other.a
			_b = other.b
			Quaternion.send(:new, @a * _a - _b.conj * @b, _b * @a + @b * _a.conj)
		elsif other.kind_of?(Numeric) && other.complex?
			# 異なるクラスでも複素数や実数なら扱える。
			# otherをQuaternionに変換してから計算してもいいが、
			# _a == other, _b == 0 であることを利用すれば効率化できる。
			Quaternion.send(:new, @a * other, @b * other.conj)
		else
			# 八元数みたいな上位の数、あるいは数でない場合を想定。
			# このクラスでは扱い方が分からないので、相手に型を揃えてもらう。
			obj_self, obj_other = other.coerce(self)
			obj_self * obj_other   # obj_selfは一般にQuaternionではなくなっている
		end
	end
end

次に演算子の左側が他クラスの場合。他クラスの定義をいじる必要は全くなく、自クラスに #coerce を実装する。戻り値の self と other の順番に注意。
例えば Complex::I * Quaternion::J という計算をすると、Complex#* は扱いの分からないクラスに困り Quaternion::J.coerce(Complex::I) と呼び出してくる。それをこちらが受けたらクラスを揃えて返すことで、結果的に相手側では Quaternion::I * Quaternion::J という計算に帰着する。

強制型変換

class Quaternion
	def coerce(other)
		if other.kind_of?(Quaternion)
			# 四元数同士。変換不要。
			[other, self]
		elsif other.kind_of?(Numeric) && other.complex?
			# 相手が複素数や実数なら、四元数に揃えることができる。
			[Quaternion.send(:new, other, 0), self]
		else
			# 相手がそれ以外なら方法が分からないので諦める。
			# ※ Rationalなどの実数クラスではsuperを呼び出す手がある。
			raise TypeError,
			      "#{other.class} can't be coerced into #{self.class}"
		end
	end
end

除算の種類

数値クラスの除算には以下の種類がある。

• #/ ： 何かしらの除算。Numeric では定義されず、各サブクラスで適切に定義する
• #div ： 余りつき除算の整数商を返す
• #fdiv ： 商（の係数）を Float で返す
• #quo ： 商をなるべく正確に返す
• Kernel.#Rational ： 整数を与えれば有理数（分数）を返すが、そうでない場合は x / y を計算して返す

除算の種類

x = Complex(4.0, 6)
y = 2
x / y          #=> (2.0+(3/1)*i)
x.div(y)       #=> NoMethodError
x.fdiv(y)      #=> (2.0+3.0i)
x.quo(y)       #=> (2.0+(3/1)*i)
Rational(x, y) #=> (2.0+(3/1)*i)

四元数の場合は、#/ は #quo のエイリアス、#div は非定義とする。#quo と #fdiv は同じ計算構造で、内部で使う除算のみが異なる。

なお Ruby 2.6 〜 2.7 にかけて、除算の結果として係数が「分母が1の Rational 」になった場合に Integer へ正規化するという処理が加わった。真面目に対応するのは大変だが、 Complex を整数の 1 で割るようにすれば同じように振る舞う、はず。

同値性

• #== ： 数値として同じならクラスが異なっていても真
• #eql? ： 相手も Quaternion で、4つの実数がそれぞれクラスを含め等しければ真
  • #hash ： オブジェクトのハッシュ値
• (Object#=== ： case式などで使う所属性判定用、デフォルトは #== と同じ)
• (BasicObject#equal? ： オブジェクトIDが等しければ真)

数値クラスで定義される同値判定には2種類ある。簡単な例では 0 と 0.0 を同じとみるかどうかが挙げられる。通常はクラス差に寛容な #== を使うはず。

#eql? はもう少し厳密で、ハッシュのキーが等しいかの判定に使われる。合わせて #hash は、#eql? が真なら同じ値を返すよう実装する必要がある。

比較

• <=> ： 大小関係を -1, 0, 1 で返す（比較不能ならnil）
• Comparable モジュールで定義済み
  • #<, #<=, #>, #>=, #between?, #clamp

複素数は大小比較できないので、これらは非定義にする…はずだった。

Ruby 2.7 で Complex#<=> が導入された²。虚部がゼロの場合は実数とみなして比較を実施するが、そうでなければ（たとえ同値でも） nil を返すという仕様になっている。実装は四則演算と同様で、唯一注意するのは相手に型を揃えてもらう前に #coerce があるか確認する（無ければ比較不能なので nil を返す）点。

インタフェース

• コンストラクタ
  • 複素数2つ・実数4つ
    • .rect (.rectangular) ： a+bj のaとbを指定する
    • .hrect (.hyperrectangular) [新規] ： w+xi+yj+zk の4つの実数を指定する
  • スカラーと3次元ベクトル
    • (名前を決めていないため未実装)
  • 極形式
    • .polar ： 極形式で指定する
  • 総合
    • Kernel.#Quaternion ： 引数に合わせて複数の形式に対応する
• 値の取得
  • 複素数2つ・実数4つ
    • #rect (#rectangular) ： 複素数の組として返す
    • #hrect (#hyperrectangular) [新規] ： 実数の組として返す
  • スカラーと3次元ベクトル
    • #real (#scalar [新規]) ： 実部を返す
    • #imag (#imaginary, #vector [新規]) ： 虚部を3次元ベクトルで返す
  • 極形式
    • #abs (#magnitude) ： 絶対値を返す
    • #arg (#angle, #phase) ： 偏角を返す
    • #axis [新規] ： 虚部の軸の方向を単位3次元ベクトルで返す
    • #polar ： 上3つをまとめて配列にして返す

Complex から何かしらの変更・拡張が必須なところ。正直なところ仕様を決め切れていない。

Complex.rect では実数2つを指定するが、四元数では実数4つ・複素数2つ・スカラー＆3次元ベクトルというパターンが考えられる。「長方形」という名前からすると対等な2つの数が良さそうなので、複素数2つを指定することにする。

Complex.polar では絶対値と偏角を指定するが、四元数では軸の情報も必要になる。極形式の指数部という類推で $\vec{n} \theta$ を指定させれば実装が楽なのだが、それでは「偏角」という名前に合わないので、偏角と軸を分離して3要素を指定することにする。

モジュール関数 Kernel.#Quaternion は、クラス内のコンストラクタならエラーを返す引数にも対応させる。具体的には、

• 文字列として与えられた引数は数値に変換する。変換は String#to_xxx より厳密で、最後まで読み取れなければエラーを返す。文法は複素数の場合が参考になる。
• できるだけ多くの数値クラスを受け入れる。例えば .rect は複素数までしか対応していないが、a+bj の式を使えば四元数や八元数に対しても妥当な結果を返せる。
  • 実際 Kernel.#Rational や Kernel.#Complex は、四元数を指定してもいいように作られている。

型変換

• Quaternion から他へ
  • #to_i, #to_r, #to_f, #to_c ： 整数・有理数・浮動小数点数・複素数に変換する
  • (#to_int ： to_i)
  • #rationalize ： #to_r と同じく有理数化するが、精度を指定可能
  • #to_q ： self を返す
  • #to_s, #inspect ： 文字列化する
• 他から Quaternion へ
  • Numeric#to_q ： 他の数値クラス（複素数までを想定）を四元数に変換する
  • String#to_q ： 文字列を解釈して四元数に変換する（極形式には未対応）
  • NilClass#to_q ： ゼロを返す
• 他から Quaternion の虚部へ
  • Numeric#j, Numeric#k [新規] ： 虚数単位を付けた数を返す

四元数を実数や複素数へ変換するときは、切り落とされる虚部が厳密にゼロ³である場合のみ許可する。

文字列化は馴染みやすさを考え "w+xi+yj+zk" という形にする（複素数の組ではない）。2つの文字列化には、各係数の変換方法および全体を括弧で囲むかの違いがある。虚部の成分が負数なら、負数の足し算ではなく正数の引き算で表すことに注意。

q = Quaternion(-1, 2, -3, 4.0) / 5
q.to_s      #=> "-1/5+2/5i-3/5j+0.8k"
q.inspect   #=> "((-1/5)+(2/5)*i-(3/5)*j+0.8k)"

String#to_q は Kernel のメソッドより緩く、文字列を解釈できるところまで読み込んで変換する。

型変換ではないが似たメソッドとして、実数を複素数の純虚数に変える Numeric#i がある。これの j, k 版があれば四元数を同様に作れるようになるので定義しておく。

虚部へ変換

1 + 2.i                 #=> (1+2i)
1 + 2.i + 3.j + 4.k     #=> (1+2i+3j+4k)
(1 + 2.i) + (3 + 4.i).j #=> (1+2i+3j+4k)

その他のメソッド

• #** ： 累乗を計算する
• #numerator, #denominator ： 四元数を「(整数係数の四元数)/(自然数)」と分数で表したときの分子・分母をそれぞれ返す

非定義にするメソッド

• 比較関連
  • #<=> （Ruby 2.6まで）
  • #<, #<=, #>, #>=, #between?, #clamp
  • #positive?, #negative?
  • #step
• 整数化
  • #ceil, #floor
  • #round
  • #truncate
• 整数商と余り
  • #div
  • #modulo (#%), #remainder
  • #divmod
• 複素数・四元数の作成
  • #i, #j [新規], #k [新規]
  • Complex#k [新規]

Numeric は実数クラスのように定義されたメソッドがあり、それらは Quaternion から呼び出すべきでない。Complex と同様に非定義にしておく。

その他

• Ruby 3.1 から、 Vector を提供する matrix ライブラリが bundled gem に降格した。そのためgemの依存関係に明示する必要が生じた。
  • 一方で Ruby 2.0 はgem版の matrix が対応していない（必須キーワード引数を使っているため）。止む無くサポートを打ち切った。
• 数年間放置していたら Complex への追従が遅れてしまったので、毎月最新のRubyを使って単体テストするようCIにスケジューリングした。

参考

• Ruby リファレンスマニュアル
  • Numeric
  • Complex
• Rubyのソースコード
  • 2.4.0: complex.c, numeric.c
  • 1.8.7: lib/complex.rb, lib/rational.rb

1. 四元数は掛ける順序が異なると結果も異なる。ここで言う除算は「除数の逆数を右から掛ける」と定義する。 ↩

2. Complex#coerce を持っているのにこれが無いことで、演算子の右辺が Complex だとエラーを起こしていた。 ↩

3. 浮動小数点数でない、ゼロ除算エラーを起こすような数。 ↩