Aizu Online Judge でC++の問題で行列の掛け算の問題があり、その時はforループを３重に使って解きました。最近、OpenGL入門という本でOpenGLの勉強をしていて、サンプルに４×４行列の掛け算をするプログラムがありました。
for (int i = 0; i < 16; ++i) {
const int j(i & 3), k(i & ~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];
}
４×４固定なので16個の配列を使います。驚いたのはiに対して3と~3のANDした結果を足して計算していること。ぱっと見てどうなっているか理解できなかったので、プログラムで一つずつ出力してみました。
まずは 3 と ~3を２進数で出力すると下記のようになります。
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で出力しました。