Numerical Recipes がダメなわけ

http://nakano.webmasters.gr.jp/nr.html


今日はNumerical Recipes(以下NP)がダメだという話を実感した。(主に僕がダメだったのだが)


数値計算の世界というのは、不具合があるかないかを判定するだけで、それなりの知識が必要とされる。
どこに境界条件があるのか、何が精度/パフォーマンスに影響を与えるのか、を判定するのはなかなか難しい。


これはユニットテストを通ればOKとかそういう世界ではなくて、整数値プログラマの世界とは違う世界だということを理解すべきなのだった。


NPの問題は

  • あんまり性能やら精度に関する考察が載ってない。
  • テーマが幅広い
  • プログラマが読むべき本に挙げられる

というところにあって、腐った実装を世に広めるのに多大に貢献しているというのがある。


せめてどれか欠けていれば問題もマシになってるはずだが、全てが揃ってしまったために、今日も「値が特殊だと精度無くなるゆっくり数値演算プログラム」が世の中の幅広い分野へはばたいていくのだった。


精度の評価ができないプログラマが数値演算プログラムを書くというのは非常に問題があって、評価方法が決まらないので最初の実装が仕様ということになってしまい、「精度が落ちるのが仕様です」となって、まともなプログラマが精度改善すると、「それバグってます」とか言われて死亡。

腐った実装を広めるのは非常に害悪、将来への負の遺産なので、それをやるくらいなら、「理論は知らないけどNetlibのライブラリ使いました」のほうがかなりマシである。まあ、つまり、NPはトータルで見てダメだと思う。



色々なテーマが載っているのはよいと思うので、この本を紹介する人は、「ただし、ここにのってるプログラムを実用するな」という言葉を添えて紹介すべきだろうと思う。

Core2は積と和ができると信じているが 

    for (i=0; i<N; i++) {
      asm volatile (""
		    "mulpd %0, %1\n\t"
		    "addpd %2, %3\n\t"
		    "mulpd %4, %5\n\t"
		    "addpd %6, %7\n\t"

		    "mulpd %0, %1\n\t"
		    "addpd %2, %3\n\t"
		    "mulpd %4, %5\n\t"
		    "addpd %6, %7\n\t"

		    "mulpd %0, %1\n\t"
		    "addpd %2, %3\n\t"
		    "mulpd %4, %5\n\t"
		    "addpd %6, %7\n\t"

		    "mulpd %0, %1\n\t"
		    "addpd %2, %3\n\t"
		    "mulpd %4, %5\n\t"
		    "addpd %6, %7\n\t"
		    "addpd %6, %7\n\t"

		    "xorpd %1, %1\n\t"
		    "xorpd %3, %3\n\t"
		    "xorpd %5, %5\n\t"
		    "xorpd %7, %7\n\t"

		    :"+x"(a),
		     "+x"(b),
		     "+x"(c),
		     "+x"(d),
		     "+x"(e),
		     "+x"(f),
		     "+x"(g),
		     "+x"(h)
		    :"m"(cache[0]));
    }

ちゃんとできる。上のが9.0 〜 9.5cycleくらい(なんかループのアラインによって変わる)
最後のaddは入れたほうが速い。なぜか。

ただ、cycleあたり2つの出力だと、4cycleでレジスタを使いきってしまい、mulpdのレイテンシ5を隠蔽できない。

ピーク性能を出すには、CoreMAのパイプラインを完全に理解し、適切な位置でロードしたりストアしたりしてあげる気遣いが大切なのである。
(実際、Core2は1cycleにひとつだけロード付き命令が実行できる(検証してない)ので不可能ではないはず)
上のはxorして日和ってる。