C++で四元数

四元数(quaternion)なるものをご存知でしょうか？
複素数をさらに拡張したような概念で虚数単位が3個もあります。

複素数は 1 と i の2つの単位を持ちますが、
四元数は 1, i, j, k の4つの単位を持ち、次式の関係があります。

ii = jj = kk = ijk = -1
ij = -ji = k,　jk = -kj = i,　ki = -ik =j

乗算に交換則は成り立ちません。

こんなん、大学でも習わんかったなあ。(僕がさぼってただけ？)
だいぶ前に、飲み会かなんかの折に、機械研のYくん(数学系)から教わりました。
こんなん、何に使うんや？　とそのときは思いましたが、
どうも3次元空間での回転をシンプルに表せるようです。
3次元との対応なら単位は3個ですむんじゃないの？ と思いますが、そうはいかないようで。

加減算、乗算、四元共役の定義

四元数 q = w + xi + yj + zk を
3次元ベクトル V = (x,y,z) をもちいて q = w + V と表すものとする。

q1 = w1 + V1, q2 = w2 + V2 に対して
加算：q1 + q2 = (w1+w2) + (V1+V2)
減算：q1 - q2 = (w1-w2) + (V1-V2)
乗算：q1 * q2 = (w1w2 - V1・V2) + (w1V2 + w2V1 + V1×V2)

また q = w + V に対し、 ~q = w - V を共役四元数という。

C++でクラス化

四元数を C++ でクラス化してみました。↓
http://licheng.sakura.ne.jp/hatena3/quaternion.lzh
演算子 +, -, * が使えます。また、.cq()で四元共役が得られます。

3次元空間での回転について

ベクトル X を、原点を中心に単位ベクトル T を回転軸として
角度 θ だけ回転させたベクトル X' を求めるには、
X,X' を 実部 = 0 の四元数とみなし、
四元数 R = cos(θ/2) + sin(θ/2)*T とすると
X' = R X ~R
とってもシンプルですね！
上にあげた C++ のクラスを使うと以下のように書けます。

#include <stdio.h>
#include <math.h>
#include "quaternion.h"

const double PI = 3.14159265358979323846;

int main(void)
{
    // もとのベクトル
    // 【例】 X1 = (100,0,0)
    quaternion X1(0, 100,0,0);
    
    // 回転軸単位ベクトル と 回転角
    // 【例】 T = (0,1,0) (つまりy軸), θ = π/2 = 90°
    vector_3D T(0,1,0);
    double th = PI / 2.0;
    
    // 回転を表す四元数
    quaternion R;
    R.SetReal (     cos(th/2.0) );
    R.SetImage( T * sin(th/2.0) );

    // 回転の演算
    quaternion X2 = R * X1 * R.cq();
    
    // 演算結果の表示
    printf("X1 = (%10.3f,%10.3f,%10.3f)\n", X1.x, X1.y, X1.z);
    printf("T  = (%10.3f,%10.3f,%10.3f)\n", T.x,  T.y,  T.z );
    printf("θ = %.3f°\n",                 th*180/PI       );
    printf("X2 = (%10.3f,%10.3f,%10.3f)\n", X2.x, X2.y, X2.z);
    
    return 0;
}