3Dグラフィックスなどで出てくるクォータニオンについて
maximaでプログラムをかいて様子を探ってみようというのが
この項のテーマです。

以下プログラム
クォータニオン

/* まず、3次の外積を定義する */
cr(a,b):=transpose(adjoint(matrix(a,b,[1,1,1])))[3];
/* ベクトル部分を取得 */
vec(a):=[a[2],a[3],a[4]];
mul(a,b):=block([ar:a[1],u:vec(a),br:b[1],v:vec(b)],return(cons(ar*br-(u . v),ar*v+br*u+cr(u,v))));

/*	回転の軸がベクトル(a,b,c),ただしa^2+b^2+c^2=1の時
	ベクトル(a,b,c)の進む向きに向かって反時計回り(右手座標系(openglなど))に
	角s回転を表す四元数
	rはqの共役四元数
	rotate(x,s,a,b,c):=mul(r(s,a,b,c),mul(x,q(s,a,b,c));
*/
q(s,a,b,c):=[cos(s/2),a*sin(s/2),b*sin(s/2),c*sin(s/2)];
r(s,a,b,c):=[cos(s/2),-a*sin(s/2),-b*sin(s/2),-c*sin(s/2)];
define(rotate(x,s,a,b,c),ratsimp(mul(r(s,a,b,c),mul(x,q(s,a,b,c)))));
define(rotate2(x,s,a,b,c),ratsimp(mul(mul(r(s,a,b,c),x),q(s,a,b,c))));

/* 結合律は成り立っているっぽい */
rotate(x,s,a,b,c)-rotate2(x,s,a,b,c);

/* x,y,z軸での回転を表現 */
define(rotatex(x,s),ratsimp(trigreduce(trigsimp(rotate(x,s,1,0,0)))));
define(rotatey(x,s),ratsimp(trigreduce(trigsimp(rotate(x,s,0,1,0)))));
define(rotatez(x,s),ratsimp(trigreduce(trigsimp(rotate(x,s,0,0,1)))));

/* 実際に試してみる。角度は右手座標で考えるのがポイント */
/* (1,0,0)をy軸を中心に90度回転すると(0,0,1) */
float(rotate([0,1,0,0],%pi/2,0,1,0));
/* (1,0,0)をz軸を中心に90度回転すると(0,-1,0) */
float(rotate([0,1,0,0],%pi/2,0,0,1));