【Unity】シェーダーでクォータニオンを扱うときの計算量について

概要

シェーダーでシミュレーションを行う際，物体の回転を扱う方法として，回転行列とクォータニオンの 2 つの選択肢があります．
回転にクォータニオンを使うと，回転行列よりも必要な変数サイズが少なくて済みますが，計算量はどうなるのか気になったので調べてみました．
シェーダーの計算負荷はアセンブリの行数と一時変数の使用量によるところが大きいらしいので，ここではアセンブリの実行行数を計算量として見ていきます．

クォータニオンによるベクトルの回転

計算式の導出

本記事ではクォータニオンを以下のように定義します．

\def\bm#1{{\boldsymbol #1}}
\begin{align}
q &= \bm{q} + q_w \\
  &= (q_x, q_y, q_z) + q_w \\
  &= (q_x, q_y, q_z, q_w)
\end{align}

また，回転を表すクォータニオンだけを扱うので，長さは１とします．

\def\bm#1{{\boldsymbol #1}}
\begin{align}
{\|q\|}^2 &= {q_x}^2 + {q_y}^2 + {q_z}^2 + {q_w}^2 = 1
\end{align}

クォータニオンによるベクトルの回転は以下のように計算されます．実行時の計算量が少なくなるように式変形します．$q^*$ は共役です．

\def\bm#1{{\boldsymbol #1}}
\begin{align}

\bm{v}'
&= q\bm{v}q^* \\
&= (\bm{q} + q_w)\bm{v}(-\bm{q} + q_w) \\
&= (\bm{q}\bm{v} + q_w\bm{v})(-\bm{q} + q_w) \\
&= -\bm{q}\bm{v}\bm{q} - q_w\bm{v}\bm{q} + q_w\bm{q}\bm{v} + {q_w}^2\bm{v} \\
&= -(-\bm{q}\cdot\bm{v} + \bm{q}\times\bm{v})\bm{q} + q_w(\bm{v}\cdot\bm{q} - \bm{v}\times\bm{q} - \bm{q}\cdot\bm{v} + \bm{q}\times\bm{v}) + {q_w}^2\bm{v} \\
&= (\bm{q}\cdot\bm{v})\bm{q} - (\bm{q}\times\bm{v})\bm{q} + 2q_w(\bm{q}\times\bm{v}) + {q_w}^2\bm{v} \\
&= (\bm{q}\cdot\bm{v})\bm{q} + (\bm{q}\times\bm{v})\cdot\bm{q} - (\bm{q}\times\bm{v})\times\bm{q} + 2q_w(\bm{q}\times\bm{v}) + {q_w}^2\bm{v} \\
&= (\bm{q}\cdot\bm{v})\bm{q} + (\bm{q}\times\bm{v})\cdot\bm{q} + (\bm{q}\cdot\bm{v})\bm{q} - (\bm{q}\cdot\bm{q})\bm{v} + 2q_w(\bm{q}\times\bm{v}) + {q_w}^2\bm{v} \\
&= 2(\bm{q}\cdot\bm{v})\bm{q} - (\bm{q}\cdot\bm{q})\bm{v} + 2q_w(\bm{q}\times\bm{v}) + {q_w}^2\bm{v} \\
&= 2((\bm{q}\cdot\bm{v})\bm{q} - (\bm{q}\cdot\bm{q})\bm{v} + q_w(\bm{q}\times\bm{v})) + (\bm{q}\cdot\bm{q} + {q_w}^2)\bm{v} \\
&= 2(\bm{q}\times(\bm{q}\times\bm{v}) + q_w(\bm{q}\times\bm{v})) + \bm{v}


\end{align}

実装

上の式を HLSL で書くとこんな感じです．

inline float3 qrot(float4 q, float3 v)
{
	float3 crs = cross(q.xyz, v);
	return v + 2 * (cross(q.xyz, crs) + q.w * crs);
}

こんなシェーダーを書いて回転行列による回転の計算量と比較します．

Shader "Test/QuaternionTest"
{
    SubShader
    {
        Pass
        {
            CGPROGRAM
            #pragma vertex vert_img
            #pragma fragment frag
            #include "UnityCG.cginc"

            float4 _Quaternion0;
            float3x3 _Rotation0;
            
            inline float3 qrot(float4 q, float3 v)
            {
                float3 crs = cross(q.xyz, v);
                return v + 2 * (cross(q.xyz, crs) + q.w * crs);
            }
            
            float4 frag (v2f_img i) : SV_Target
            {
                float3 p = i.pos.xyz;
            #if 1
                p = qrot(_Quaternion0, p);
            #else
                p = mul(_Rotation0, p);
            #endif
                return float4(p, 1);
            }
            ENDCG
        }
    }
}

コンパイルしてアセンブリの該当部分を見てみます．

クォータニオンの場合の計算量

アセンブリの該当部分は以下です．

0: mul r0.xyz, v0.yzxy, cb0[2].zxyz
1: mad r0.xyz, cb0[2].yzxy, v0.zxyz, -r0.xyzx
2: mul r1.xyz, r0.yzxy, cb0[2].zxyz
3: mad r1.xyz, cb0[2].yzxy, r0.zxyz, -r1.xyzx
4: mad r0.xyz, cb0[2].wwww, r0.xyzx, r1.xyzx
5: mad o0.xyz, r0.xyzx, l(2.000000, 2.000000, 2.000000, 0.000000), v0.xyzx
6: mov o0.w, l(1.000000)
7: ret 

6 ステップかかっています．0, 1 が cross(q.xyz, v) 部分，2, 3 が cross(q.xyz, crs) 部分，4 が (... + q.w * crs) 部分，5 が v + 2 * (...) 部分で，6 は float4(p, 1) 部分です．
レジスタは v0 と o0 を除いて 3 つ必要なようです．クォータニオン自身 cb0[2] で 1 つ分，一時変数 r0 r1 で 2 つ分です．

回転行列の場合の計算量

アセンブリの該当部分は以下です．

0: mul r0.xyz, v0.yyyy, cb0[4].xyzx
1: mad r0.xyz, cb0[3].xyzx, v0.xxxx, r0.xyzx
2: mad o0.xyz, cb0[5].xyzx, v0.zzzz, r0.xyzx
3: mov o0.w, l(1.000000)
4: ret 

3 ステップかかっています．純粋に行列とベクトルの積の計算をしていますね．
レジスタは v0 と o0 を除いて 4 つ必要なようです．回転行列自身 cb0[3] cb0[4] cb0[5] で 3 つ分，一時変数 r0 で 1 つ分です．

考察

クォータニオンによる回転は，回転行列に比べて，計算量は 2 倍，レジスタの使用数は 0.75 倍，のようなので，よほど多種な回転を同時に保持するような状況でない限り，回転行列を使ったほうがよさそうです．
ただし，シミュレーションで回転を積分する場合は適度に normalize() する必要があります．クォータニオンは 1 回で済むのに対し，回転行列は各軸それぞれに正規化が必要なので 3 回 normalize() が必要になります．normalize() の計算量は 3 ステップなので，クォータニオンの正規化は 3 ステップ，回転行列の正規化は 9 ステップの計算量になります．

float4 frag (v2f_img i) : SV_Target
{
    float3 p = i.pos.xyz;
    p = normalize(p);
    return float4(p, 1);
}

0: dp3 r0.x, v0.xyzx, v0.xyzx
1: rsq r0.x, r0.x
2: mul o0.xyz, r0.xxxx, v0.xyzx
3: mov o0.w, l(1.000000)
4: ret 

そのため，計算テクスチャから回転を取り出したあと，正規化してからベクトルの回転を行う場合は，クォータニオンは 9 ステップ，回転行列は 12 ステップ，それぞれかかることになります．

float4 frag (v2f_img i) : SV_Target
{
    float3 p = i.pos.xyz;
#if 1
    float4 q = normalize(_Quaternion0);
    p = qrot(q, p);
#else
    float3x3 R;
    R._m00_m10_m20 = normalize(_Rotation0._m00_m10_m20);
    R._m01_m11_m21 = normalize(_Rotation0._m01_m11_m21);
    R._m02_m12_m22 = normalize(_Rotation0._m02_m12_m22);
    p = mul(R, p);
#endif
    return float4(p, 1);
}

クォータニオンから回転行列への変換

回転行列のほうが回転自体の計算量は少ないですが，回転行列（9 要素）はクォータニオン（4 要素）に比べてサイズが大きいので，計算テクスチャで回転を保存する場合はクォータニオンを使いたいです．
計算テクスチャからクォータニオンを取り出し，回転行列に変換してから回転させることを考えます．

クォータニオンから回転行列への変換

変換式は以下です．

\def\bm#1{{\boldsymbol #1}}
R = 
\left(
  \begin{array}
    1-2({q_y}^2+{q_z}^2) & 2(q_xq_y-q_wq_z) & 2(q_xq_z+q_wq_y) \\\
    2(q_xq_y+q_wq_z) & 1-2({q_x}^2+{q_z}^2) & 2(q_yq_z-q_wq_x) \\\
    2(q_xq_z-q_wq_y) & 2(q_yq_z+q_wq_x) & 1-2({q_x}^2+{q_y}^2) \\
  \end{array}
\right)

ただ，これだと計算量が増えるので，導出の途中から派生します．

$R$ の軸 $\boldsymbol{e}_x',\boldsymbol{e}_y',\boldsymbol{e}_z'$ を求めます．冒頭で導出した式を使います．

\def\bm#1{{\boldsymbol #1}}
\begin{align}

\bm{e}_x'
&= q^*\bm{e}_xq \\
&= \bm{e}_x + 2(\bm{q}\times(\bm{q}\times\bm{e}_x) + q_w(\bm{q}\times\bm{e}_x)) \\

\bm{e}_y'
&= q^*\bm{e}_yq \\
&= \bm{e}_y + 2(\bm{q}\times(\bm{q}\times\bm{e}_y) + q_w(\bm{q}\times\bm{e}_y)) \\

\bm{e}_z'
&= \bm{e}_x' \times \bm{e}_y' \\

\end{align}

ここから効率的な計算方法を考えます．$\boldsymbol{e}_x,\boldsymbol{e}_y$ はリテラルなので展開しておきます．

\def\bm#1{{\boldsymbol #1}}
\begin{align}
\bm{q}\times\bm{e}_x &= (q_x,q_y,q_z) \times (1,0,0) = (0,q_z,-q_y) \\
\bm{q}\times\bm{e}_y &= (q_x,q_y,q_z) \times (0,1,0) = (-q_z,0,q_x) \\
\end{align}

float4 のベクトル $\boldsymbol{f}$ を以下のように定義します．

\def\bm#1{{\boldsymbol #1}}
\bm{f} = (-2q_x, -2q_y, 2q_z, 0)

$\boldsymbol{q}\times\boldsymbol{e}_x$ と $\boldsymbol{q}\times\boldsymbol{e}_x$ を $\boldsymbol{f}$ で書き直します．

\def\bm#1{{\boldsymbol #1}}
\begin{align}
\bm{q}\times\bm{e}_x &= \frac{1}{2}\bm{f}_{wzy} \\
\bm{q}\times\bm{e}_y &= -\frac{1}{2}\bm{f}_{zwx} \\
\end{align}

$R$ の軸を書き直します．

\def\bm#1{{\boldsymbol #1}}
\begin{align}

\bm{e}_x' &= \bm{q}\times\bm{f}_{wzy} + q_w\bm{f}_{wzy} + \bm{e}_x \\
\bm{e}_x' &= -\bm{q}\times\bm{f}_{zwx} - q_w\bm{f}_{zwx} + \bm{e}_y \\
\bm{e}_z' &= \bm{e}_x' \times \bm{e}_y' \\

\end{align}

これで $R$ が求まりました．

実装

以上を HLSL で実装します．極力計算量が減るように cross() も展開しています．可読性は最悪ですが効率を優先します．

inline float3x3 quaternionToMatrix(float4 q)
{
    float4 f = q * float4(-2, -2, 2, 0);
    float3 x = float3(1, 0, 0) + q.yzx * f.ywz + q.zxy * f.zyw - q.w * f.wzy;
    float3 y = float3(0, 1, 0) - q.yzx * f.xzw - q.zxy * f.wxz + q.w * f.zwx;
    float3 z = cross(x, y);

    float3x3 R;
    R._m00_m10_m20 = x;
    R._m01_m11_m21 = y;
    R._m02_m12_m22 = z;
    
    return R;
}

シェーダーに組み込みます．

Shader "Test/QuaternionTest"
{
    SubShader
    {
        Pass
        {
            CGPROGRAM
            #pragma vertex vert_img
            #pragma fragment frag
            #include "UnityCG.cginc"

            float4 _Quaternion0;
            
            inline float3x3 quaternionToMatrix(float4 q)
            {
                float4 f = q * float4(-2, -2, 2, 0);
                float3 x = float3(1, 0, 0) + q.yzx * f.ywz + q.zxy * f.zyw - q.w * f.wzy;
                float3 y = float3(0, 1, 0) - q.yzx * f.xzw - q.zxy * f.wxz + q.w * f.zwx;
                float3 z = cross(x, y);

                float3x3 R;
                R._m00_m10_m20 = x;
                R._m01_m11_m21 = y;
                R._m02_m12_m22 = z;
                
                return R;
            }
            
            float4 frag (v2f_img i) : SV_Target
            {
                float3 p = i.pos.xyz;

                float3x3 rot = quaternionToMatrix(_Quaternion0);
                p = mul(rot, p);

                return float4(p, 1);
            }
            ENDCG
        }
    }
}

コンパイルしてアセンブリを見てみます．

   0: mul r0.xyzw, cb0[2].xyzw, l(-2.000000, -2.000000, 2.000000, 0.000000)
   1: mad r1.xyz, cb0[2].yzxy, r0.ywzy, l(1.000000, 0.000000, 0.000000, 0.000000)
   2: mad r1.xyz, cb0[2].zxyz, r0.zywz, r1.xyzx
   3: mad r1.xyz, -cb0[2].wwww, r0.wzyw, r1.xyzx
   4: mad r2.xyz, -cb0[2].yzxy, r0.xzwx, l(0.000000, 1.000000, 0.000000, 0.000000)
   5: mad r2.xyz, -cb0[2].zxyz, r0.wxzw, r2.xyzx
   6: mad r0.xyz, cb0[2].wwww, r0.zwxz, r2.xyzx
   7: mul r2.xyz, r0.yzxy, r1.zxyz
   8: mad r2.xyz, r1.yzxy, r0.zxyz, -r2.xyzx
   9: mul r0.xyz, r0.xyzx, v0.yyyy
  10: mad r0.xyz, r1.xyzx, v0.xxxx, r0.xyzx
  11: mad o0.xyz, r2.xyzx, v0.zzzz, r0.xyzx
  12: mov o0.w, l(1.000000)
  13: ret 

9 以降が行列とベクトルの積なので，0 ~ 8 までの 9 ステップが回転行列の計算のようです．それぞれ，0 が f = ... の行，1 ~ 3 が x = ... の行，4 ~ 6 が y = ... の行，7 ~ 8 が z = ... の行です．

考察

クォータニオンでの回転が 6 ステップ，行列での回転が 3 ステップ，クォータニオンから回転行列への変換が 9 ステップなので，同一の回転を使って 4 回以上回転を行う場合は回転行列に変換したほうがよさそうです．

クォータニオンによる回転の回転

クォータニオンの回転はクォータニオン同士の積です．

\def\bm#1{{\boldsymbol #1}}
\begin{align}
q_0q_1 
&= (\bm{q}_0 + {q_0}_w)(\bm{q}_1 + {q_1}_w) \\
&= - \bm{q}_0\cdot\bm{q}_1 + \bm{q}_0\times\bm{q}_1 + {q_1}_w\bm{q}_0 + {q_0}_w\bm{q}_1 + {q_0}_w{q_1}_w \\
&= (\bm{q}_0\times\bm{q}_1 + {q_1}_w\bm{q}_0 + {q_0}_w\bm{q}_1) + ({q_0}_w{q_1}_w -\bm{q}_0\cdot\bm{q}_1)
\end{align}

実装

HLSL で実装します．アセンブリでなるべく mad が使われるよう足し算の順番に気を付けます．

inline float4 qmul(float4 a, float4 b)
{
    return float4((cross(a.xyz, b.xyz) + a.w * b.xyz) + b.w * a.xyz, a.w * b.w - dot(a.xyz, b.xyz));
}

シェーダーに組み込みます．アセンブリで「回転の回転」の部分が見やすいように，該当部分を if 節にしています．

Shader "Test/QuaternionTest"
{
    SubShader
    {
        Pass
        {
            CGPROGRAM
            #pragma vertex vert_img
            #pragma fragment frag
            #include "UnityCG.cginc"

            float4 _Quaternion0;
            float4 _Quaternion1;
            float3x3 _Rotation0;
            float3x3 _Rotation1;
            
            inline float4 qmul(float4 a, float4 b)
            {
                return float4((cross(a.xyz, b.xyz) + a.w * b.xyz) + b.w * a.xyz, a.w * b.w - dot(a.xyz, b.xyz));
            }
            
            inline float3 qrot(float4 q, float3 v)
            {
                float3 crs = cross(q.xyz, v);
                return v + 2 * (cross(q.xyz, crs) + q.w * crs);
            }
            
            float4 frag (v2f_img i) : SV_Target
            {
                float3 p = i.pos.xyz;
            #if 1
                float4 q = float4(0, 0, 0, 1);
                if (i.uv.x > 0.5)
                {
                    q = qmul(_Quaternion0, _Quaternion1);
                }
                p = qrot(q, p);
            #else
                float3x3 R = float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
                if (i.uv.x > 0.5)
                {
                    R = mul(_Rotation0, _Rotation1);
                }
                p = mul(R, p);
            #endif
                return float4(p, 1);
            }
            ENDCG
        }
    }
}

コンパイルしてアセンブリを見てみます．

クォータニオンの場合の計算量

アセンブリの該当部分は以下です．

 0: mul r0.xyz, cb0[2].yzxy, cb0[3].zxyz
 1: mad r0.xyz, -cb0[2].zxyz, cb0[3].yzxy, r0.xyzx
 2: mad r0.xyz, cb0[2].wwww, cb0[3].xyzx, r0.xyzx
 3: mad r0.xyz, cb0[3].wwww, cb0[2].xyzx, r0.xyzx
 4: dp3 r1.x, cb0[2].xyzx, cb0[3].xyzx
 5: mad r0.w, cb0[2].w, cb0[3].w, -r1.x
 6: lt r1.x, l(0.500000), v1.x
 7: movc r0.xyzw, r1.xxxx, r0.xyzw, l(0,0,0,1.000000)
 8: mul r1.xyz, r0.zxyz, v0.yzxy
 9: mad r1.xyz, r0.yzxy, v0.zxyz, -r1.xyzx
10: mul r2.xyz, r0.zxyz, r1.yzxy
11: mad r0.xyz, r0.yzxy, r1.zxyz, -r2.xyzx
12: mad r0.xyz, r0.wwww, r1.xyzx, r0.xyzx
13: mad o0.xyz, r0.xyzx, l(2.000000, 2.000000, 2.000000, 0.000000), v0.xyzx
14: mov o0.w, l(1.000000)
15: ret 

0 ~ 5 が qmul() の計算，6 ~ 7 が if () の部分（q に初期値を入れるか qmul() の結果を入れるか選択している）で，8 ~ 13 が qrot() の計算，なので，クォータニオン積の計算には 6 ステップかかっています．

回転行列の場合の計算量

アセンブリの該当部分は以下です．if 節のせいでややこしくなっています．

 0: mov o0.w, l(1.000000)
 1: mul r0.xyz, cb0[5].xyzx, cb0[7].yyyy
 2: mad r0.xyz, cb0[4].xyzx, cb0[7].xxxx, r0.xyzx
 3: mad r0.xyz, cb0[6].yxzy, cb0[7].zzzz, r0.yxzy
 4: mov r1.x, r0.y
 5: mul r2.xyz, cb0[5].xyzx, cb0[8].yyyy
 6: mad r2.xyz, cb0[4].xyzx, cb0[8].xxxx, r2.xyzx
 7: mad r2.xyz, cb0[6].xzyx, cb0[8].zzzz, r2.xzyx
 8: mov r1.y, r2.x
 9: mul r3.xyz, cb0[5].xyzx, cb0[9].yyyy
10: mad r3.xyz, cb0[4].xyzx, cb0[9].xxxx, r3.xyzx
11: mad r3.xyz, cb0[6].xyzx, cb0[9].zzzz, r3.xyzx
12: mov r1.z, r3.x
13: lt r0.w, l(0.500000), v1.x
14: movc r1.xyz, r0.wwww, r1.xyzx, l(1.000000,0,0,0)
15: dp3 o0.x, r1.xyzx, v0.xyzx
16: mov r2.x, r0.z
17: mov r0.y, r2.z
18: mov r0.z, r3.y
19: mov r2.z, r3.z
20: movc r1.xyz, r0.wwww, r2.xyzx, l(0,0,1.000000,0)
21: movc r0.xyz, r0.wwww, r0.xyzx, l(0,1.000000,0,0)
22: dp3 o0.y, r0.xyzx, v0.xyzx
23: dp3 o0.z, r1.xyzx, v0.xyzx
24: ret 

0 は return float4(p, 1) の 1 の部分です．
1 ~ 3, 5 ~ 7, 9 ~ 11 が行列積 mul() の計算です．
13 ~ 14, 20 ~ 21 が if () の部分です．R に初期値を入れるか mul() の結果を入れるか判断しています．
4, 8, 12, 16 ~ 19 は行列の転置を取っています．
15, 22 ~ 23 が行列とベクトルの積 mul(R, p) の計算です．

純粋に回転行列の行列積にかかる計算量は 9 ステップですが，行列の取り回しにも計算量が割かれることを考えると，実質的な計算量はもっとかかりそうです．
ちなみに，今回の場合だと，mul(R, p) の計算方法をべた書きしておくと転置の計算はなくなります．汎用性はないですが．．．

// p = mul(R, p);
p = R._m00_m10_m20 * p.x + R._m01_m11_m21 * p.y + R._m02_m12_m22 * p.z;

 0: mul r0.xyz, cb0[5].xyzx, cb0[8].yyyy
 1: mad r0.xyz, cb0[4].xyzx, cb0[8].xxxx, r0.xyzx
 2: mad r0.xyz, cb0[6].xyzx, cb0[8].zzzz, r0.xyzx
 3: lt r0.w, l(0.500000), v1.x
 4: movc r0.xyz, r0.wwww, r0.xyzx, l(0,1.000000,0,0)
 5: mul r0.xyz, r0.xyzx, v0.yyyy
 6: mul r1.xyz, cb0[5].xyzx, cb0[7].yyyy
 7: mad r1.xyz, cb0[4].xyzx, cb0[7].xxxx, r1.xyzx
 8: mad r1.xyz, cb0[6].xyzx, cb0[7].zzzz, r1.xyzx
 9: movc r1.xyz, r0.wwww, r1.xyzx, l(1.000000,0,0,0)
10: mad r0.xyz, r1.xyzx, v0.xxxx, r0.xyzx
11: mul r1.xyz, cb0[5].xyzx, cb0[9].yyyy
12: mad r1.xyz, cb0[4].xyzx, cb0[9].xxxx, r1.xyzx
13: mad r1.xyz, cb0[6].xyzx, cb0[9].zzzz, r1.xyzx
14: movc r1.xyz, r0.wwww, r1.xyzx, l(0,0,1.000000,0)
15: mad o0.xyz, r1.xyzx, v0.zzzz, r0.xyzx
16: mov o0.w, l(1.000000)
17: ret 

考察

クォータニオン積の計算が 6 ステップ，行列積の計算が 9 ステップ，なので，回転の回転を計算したい場合はクォータニオンが有利です．また，行列の取り回しも考えると，回転行列を使う場合はもっと計算量が多くなりそうです．
しかし，ベクトルを２回回転したいだけなら，クォータニオンの場合は 6 x 2 = 12 ステップ，回転行列の場合は 3 x 2 = 6 ステップなので，回転行列のほうが有利です．
どちらのほうが計算量を少なくできるかは目的によって変わってきそうです．

まとめ

計算量をまとめます．

計算	式	ステップ数
クォータニオンでのベクトルの回転	$q\boldsymbol{v}q^*$	6
回転行列でのベクトルの回転	$R\boldsymbol{v}$	3
クォータニオンでのクォータニオンの回転	$q_0q_1$	6
回転行列での回転行列の回転	$R_0R_1$	9
クォータニオンの正規化	$normalize(q)$	3
回転行列の正規化	$normalize(R)$	9
クォータニオンから回転行列への変換	$R(q)$	9

純粋にベクトルの回転であれば回転行列が最強です．それ以外の操作が多いならクォータニオンのほうがよさそうです．

雑記

クォータニオンの計算は，最初はネットで拾ってきた関数をそのまま使っていましたが，計算手順を工夫したら計算量が半減して驚きました．シェーダーでもリファクタリングは大事ですね．

クォータニオンの計算はこのあたりが詳しいです．