SSE


Twitterで渋川さんに煽られて、そういえばそういうのは全く興味がなく知らなかったので調べてみた。libpngにコミットするとか実に気の遠い話であるがとりあえず手がかりに調べたりしてみた。もともとIntelx86の独自拡張命令の一環で浮動小数点の計算をちょっと速くするシロモノだ、くらいにしか思っていなかったのだが、要は128ビットのレジスタを複数用意してベクトル計算のようなものを実行できるらしい。スカラーvsベクトルCPUの宗教論争がいかに無意味で、量産によるコスト低減が最も大切だった…的な。

基本的には、浮動小数点を計算するための128ビットレジスタなので、floatが4コ並んでいるか、doubleが2コ並んでいるかどっちかで使えばよい。前者の場合は __m128 という型のついた変数を宣言して使えばよいし、後者の場合は __m128d という型をつけて使えばよい。 __m128i というのあって、これは128bit整数というよりは、128ビットのフラグの並びとかそういう考えをすればよい…のかな。基本的な関数は _mm_spam_pd を使う。load, storeで値を入れる。あとはadd, mul, div, subとか。cmpもいくつかあるがヤヤコシイな…

で、ネイピア数を適当に求める関数

double f0(unsigned long long n){
  double sum = 1.0;
  double fact = 1.0;
  unsigned long long i;
  for(i=2; i<=n; ++i){
    sum += 1.0 / fact;
    fact *= i;
  }
  return sum;
}

と、SSEを使って計算する関数

double f1(unsigned long long n){
  double d[2] = {1.0, 1.0}; // even, odd                                                                                                        
  double f[2] = {2.0, 6.0}; // 2n!, (2n+1)!                                                                                                     
  double ns[2] = {3.0, 4.0}; // n, n+1                                                                                                          
  double one[2] = {1.0, 1.0}; // 1, 1                                                                                                           
  __m128d sums = _mm_load_pd(d);
  __m128d fact = _mm_load_pd(f);
  __m128d is = _mm_load_pd(ns);
  __m128d ones = _mm_load_pd(one);
  unsigned long long i;
  for(i=1; i<=n; i+=2){
    sums = _mm_add_pd(sums, _mm_div_pd(ones, fact));
    fact = _mm_mul_pd(fact, is);
    is = _mm_add_pd(is, ones);
    fact = _mm_mul_pd(fact, is);
    is = _mm_add_pd(is, ones);
  }
  _mm_store_pd(d, sums);
  return d[0] + d[1];
}

の性能を比べてみたら

$ cc napier.c -mfpmath=sse -msse4 && ./a.out                                                                                                  
(0 msec) const:         2.71828182845904509080                                                                                                
(1894.87msec) f0:       2.71828182845904553488                                                                                                
(1224.71msec) f1:       2.71828182845904553488  

となって、まあ2/3くらいの性能差が出た*1。ベクトルっぽく計算を書くのに割と工夫が必要で面倒臭かったわりに1.5倍くらいしか速くならないので使いドコロは難しいなあと思いつつも、世界中のコンピュータで動くなら…と思ったり思わなかったり。
ひとつわからないのはgccのドキュメントを見る限り命名規則が全く違っていて__builtin_i32_hamみたいなのが多い。この _mm_spam_pd ってのはどこからくるのかちょっと分からないで

どうでもいいけどARM NEONはどうしてこんなに拡張命令が充実してるんだろう。やっぱそういうのがARMがATOMに勝ちそうな理由なのかなあ。

*1:cpuinfoには "Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz"とか書いてあった