4×4行列の掛け算のプログラム

Aizu Online Judge でC++の問題で行列の掛け算の問題があり、その時はforループを3重に使って解きました。最近、OpenGL入門という本でOpenGLの勉強をしていて、サンプルに4×4行列の掛け算をするプログラムがありました。

for (int i = 0; i < 16; ++i) {
    const int j(i &amp; 3), k(i &amp; ~3);
    t[i] =
       a[ 0 + j] * b[k + 0] +
       a[ 4 + j] * b[k + 1] +
       a[ 8 + j] * b[k + 2] +
       a[12 + j] * b[k + 3];
}

4×4固定なので16個の配列を使います。驚いたのはiに対して3と~3のANDした結果を足して計算していること。ぱっと見てどうなっているか理解できなかったので、プログラムで一つずつ出力してみました。

まずは 3 と ~3を2進数で出力すると下記のようになります。

 3 : 00000000000000000000000000000011
~3 : 11111111111111111111111111111100

これをループで0 ~ 15を出力すると、& 3は

 0 &  3 = 0
 1 &  3 = 1
 2 &  3 = 2
 3 &  3 = 3
 4 &  3 = 0
 5 &  3 = 1
 6 &  3 = 2
 7 &  3 = 3
 8 &  3 = 0
 9 &  3 = 1
10 &  3 = 2
11 &  3 = 3
12 &  3 = 0
13 &  3 = 1
14 &  3 = 2
15 &  3 = 3

& ~3は

 0 & ~3 = 0
 1 & ~3 = 0
 2 & ~3 = 0
 3 & ~3 = 0
 4 & ~3 = 4
 5 & ~3 = 4
 6 & ~3 = 4
 7 & ~3 = 4
 8 & ~3 = 8
 9 & ~3 = 8
10 & ~3 = 8
11 & ~3 = 8
12 & ~3 = 12
13 & ~3 = 12
14 & ~3 = 12
15 & ~3 = 12

となります。0〜15は配列の下記のように割り当てます。

 0  1  2  3 
 4  5  6  7 
 8  9 10 11 
12 13 14 15

行列の掛け算は横✕縦でかけていきますので、式を出力するとこのようになります。

C[ 0]=A[ 0]*B[ 0]+A[ 4]*B[ 1]+A[ 8]*B[ 2]+A[12]*B[ 3]
C[ 1]=A[ 1]*B[ 0]+A[ 5]*B[ 1]+A[ 9]*B[ 2]+A[13]*B[ 3]
C[ 2]=A[ 2]*B[ 0]+A[ 6]*B[ 1]+A[10]*B[ 2]+A[14]*B[ 3]
C[ 3]=A[ 3]*B[ 0]+A[ 7]*B[ 1]+A[11]*B[ 2]+A[15]*B[ 3]
C[ 4]=A[ 0]*B[ 4]+A[ 4]*B[ 5]+A[ 8]*B[ 6]+A[12]*B[ 7]
C[ 5]=A[ 1]*B[ 4]+A[ 5]*B[ 5]+A[ 9]*B[ 6]+A[13]*B[ 7]
C[ 6]=A[ 2]*B[ 4]+A[ 6]*B[ 5]+A[10]*B[ 6]+A[14]*B[ 7]
C[ 7]=A[ 3]*B[ 4]+A[ 7]*B[ 5]+A[11]*B[ 6]+A[15]*B[ 7]
C[ 8]=A[ 0]*B[ 8]+A[ 4]*B[ 9]+A[ 8]*B[10]+A[12]*B[11]
C[ 9]=A[ 1]*B[ 8]+A[ 5]*B[ 9]+A[ 9]*B[10]+A[13]*B[11]
C[10]=A[ 2]*B[ 8]+A[ 6]*B[ 9]+A[10]*B[10]+A[14]*B[11]
C[11]=A[ 3]*B[ 8]+A[ 7]*B[ 9]+A[11]*B[10]+A[15]*B[11]
C[12]=A[ 0]*B[12]+A[ 4]*B[13]+A[ 8]*B[14]+A[12]*B[15]
C[13]=A[ 1]*B[12]+A[ 5]*B[13]+A[ 9]*B[14]+A[13]*B[15]
C[14]=A[ 2]*B[12]+A[ 6]*B[13]+A[10]*B[14]+A[14]*B[15]
C[15]=A[ 3]*B[12]+A[ 7]*B[13]+A[11]*B[14]+A[15]*B[15]

検証に使ったプログラムはこれです。clang++ std-+11で出力しました。