Even命令が大幅に過剰

みんなそうだと思うけど、even命令が大幅に過剰になっちまった。
しかも、僕のコードは1ワード1ロードシャッフル無しだから、それはもう多すぎて

http://d.hatena.ne.jp/kikx/20090120

これはループのアンロール前で内側に14回のループが残ってたんだけど、
こんなんで44倍だったからアンロールしてoddに移動できるのを順番に移動してったんですよ。

(y >> 1) ^ mag01[y & 1]
y << 7
y << 15

これらは上から、

even*1 + odd*3
odd*1
odd*2

を全部やってもまだ余ってて、しかたがないから

y >> 11

の半分くらいをodd*3で置き換えて、69倍速になったのです。

せっかくのgather命令

バイト単位じゃないともったいない。

z = si_lqx(spu_slqw(spu_gather(y), 4), mag_lut);
r = spu_xor(spu_rlmaskqw(y,-1), z);

これだと、ワード単位のgatherになっちゃって結果が4ビットしか得られないんですよ。
1命令使って、結果がたったの4ビットですよ。
シャッフルしてからバイト単位gatherすれば、結果は16ビットになってうはうはだと思ったんです。

x = spu_shuffle(y, _, s);
z = si_lqx(spu_gather((vector byte)x), mag_lut);
r = spu_xor(spu_rlmaskqw(y,-1), z);

シャッフルで下4バイトに0を放り込んで、その上に目的のバイトを並べとけば4ビットシフトはいらなくなりました。
しかも、シャッフルの引数が一個余ってる。

シャッフルの引数を有効利用するためにふたつを一度にやってみると

r0 = (y0 >> 1) ^ mag01[y0 & 1]
r1 = (y1 >> 1) ^ mag01[y1 & 1]

は

x = spu_gather((vector byte)spu_shuffle(y0, y1, s));
z0 = si_lqx(x, mag_lut0);
z1 = si_lqx(x, mag_lut1);
r0 = spu_xor(spu_rlmaskqw(y0,-1), z0); 
r1 = spu_xor(spu_rlmaskqw(y1,-1), z1); 

これで、even*2+odd*6だからひとつ当たりだとeven*1+odd*3になりました。
ルックアップテーブルが大きくなりますが、256エントリのテーブルをふたつは余裕で収まったので問題なし。

ビット転置

69倍になったあとは、一ヶ月ずっとビット転置のやり方を考えてました。
ワード単位のときのような頭のおかしい配置方法があるんじゃないかと探してたんです。

結局みつからなかったので、頭から順番にビットを詰めていく方法になりました。

qword mt_bs[32][5];

の mt_bs[i][0] には 0から127番目のワードのiビット目を並べます。残りも同じように詰めます。
実際に書いてみると分かるのですが

y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];

このkk+1が鬼門で、kk+1ワード目にあったビットをkkワード目の位置までもってくるのに、1ビットシフトが必要なだけでなく、足りない分を隣のqwordから1ビットもらって来ないといけません。
31ビット分同じことをやらないといけないので、31回のrlqwとselbが必要になります。

ここで、MTの状態空間は本当に624ワードだったかを思い出してみます。
松本先生の書いたドキュメントを読んでると、状態空間は624ワードじゃなくて623ワード+1ビットであることが分かります。
これをワード単位に切り上げた実装が624ワードになってるだけです。
623ワード+1ビットで更新処理を書いてみると、

word mt[623];
word last_bit;

for (kk = 0; kk < 623; kk++) {
  y = last_bit | (mt[kk] & LOWER_MASK);
  last_bit = mt[kk] & UPPER_MASK;
  mt[kk] = mt[(kk+M-1) % 623] ^ (y >> 1) ^ mag01[y & 0x1UL];
}

これを元に実装すると、1ビットシフトは31回も必要なくて、最上位ビット1回だけになります。