a*b/65536(ただしa=[-255, 255], b=[0, 65535])をpmulhwで計算するとき、pmulhwは符号付き16bit乗算命令だから工夫が必要だ。aは符号付き16bitに収まるから問題ないとして、bは[32768, 65535]の範囲で符号付き16bitに収まらない。bは符号「なし」16bitには収まるので、符号付き16bit×符号なし16bitの乗算を行えばよいことになる。でもpmulhw命令は二つの項を符号付きとして扱ってしまう。bの範囲[32768, 65535]は16進数では[8000h, ffffh]となり、二の補数表現で符号付き16bitと解釈すると[-32768, -1]となる。例えば255×65535を計算したいとき、pmulhwはこれを255×-1として解釈してしまう。すると 255×65535 = 16711425 = 00feff01h、255×-1 = -255 = ffffff01h のように当然答えは食い違ってしまう。pmulhwは答え32bitのうち上位16bitを求める命令なので、真の値は00feh、間違った値はffffhとなる。その差はffh。
符号付き乗算は四つの場合に分けられる。aが正でbが正として扱われるとき、aが正でbが負として扱われるとき、aが負でbが正として扱われるとき、aが負でbが負として扱われるときの四つだ。aとbが正のときは問題がない。aが負でbが正のときも問題がない。問題はbが負のとき(負として扱われるとき)である。
で、結局どうすればいいのかというと、ハッカーのたのしみにも書いてあるのだけど、bが負として扱われたときに間違った値の上位16bitにaを加算してやればいい。前の例だとffhの差があったのはaが255だったからと言うわけ。
というわけで、a*b/65536を求めるコードは次のようになる(mm0=a四つ分、mm1=b四つ分、結果はmm0)。
movq mm2, mm1
psraw mm2, 15 ;mm2をbの符号で埋める。
pand mm2, mm0 ;bが負ならa、そうでないなら0。
pmulhw mm0, mm1 ;mm0=a×b (bを符号付きとして扱った場合)
paddw mm0, mm2 ;mm0=a×b (bを符号なしとして扱った場合)
とまあ、単にかけ算をしたいだけでこういう事を分かっていないといけない。もっと言えば証明できないといけない。みんなこういうのってどこで覚えるのかなぁ。