憂鬱な午後のひととき


10月30日_
眠たくて頭痛いけど, せっかく伊庭先生がMCMCについて話してくださるというのに参加しないという選択枝は無い. こないだ社内勉強会で川合さんによる統計熱力学でのモンテカルロの話を聞いたばかりだしな. というわけで会社へ.
次の魔方陣の個数を拡張アンサンブルモンテカルロ法を用いて求める話もとても面白かった. 他の話も聞きたかったけど残念ながらここで帰宅.
10月29日_
恐ろしく頭が痛い. あまりに痛くて全然寝られない. 明け方漸くうとうとした.
10月28日_

インテルソフトウェア・カンファレンス2011に行ってきた.

無料なのに軽食が出て, それがえらくおいしかった. というのはおいといて.
Cilk Plusの配列構文はおもろいな.

 C[0:n] = A[0:n] + B[0:n];
 if (A[:] > 0) { B[:] += C[:]; }

みたいなの. 関係ないけど最初"Click Plus"と読んでて「なんでシリックって発音してるのかな」と思ってたのは内緒.

parallel_forはさりげなくVC2010でもppl.hでさりげなく対応してたりするし今度ちゃんと使ってみよう.

コンパイラにさまざまな情報を与えられる#pragmaや__assume文. VCのはいまいち(書いてもコード生成に殆ど影響なさげ)だけど, iccのは相当強力らしい.

__assume_aligned(p, 32); // pが32byte alignmentされている
__assume(n % 8 == 0);
#pragma loop count avg(32); // ループの平均回数が32回

VTuneの使い方のデモもほどよくディープでよかった. float→intキャストが遅かった時代のasm { fld x; fistp ix; }はよく使ったなあ.

期待以上のカンファレンスだった.

10月27日_
朝, 急にオムライスが食べたくなって作ってみた. 初めてだが意外にうまくできておいしかった. これならまた作ってもよい.
10月25日_
『結婚小説』(中島たい子)
こないだ読んだ『漢方小説』が面白かったので他のも読んでみる. 結婚をテーマにした小説を書けと命じられたヒロインが「したこともないのに書けるか」と怒りながら体当たり取材をする. ところどころに挿入される当事者なのに第三者的視点の結婚取材メモが鋭い. ミイラ取りがミイラになるか. ありきたりなハッピーエンドではないところがいかにも著者らしい.
10月20日_
・icacheサイズとinline展開の兼ね合いが難しい. 基本的に手元のSandy Bridgeでやってるといつのまにか一世代前のXeonだと遅くなってたりして困る.
・自宅サーバのkernelを3.0.4にアップグレード. もともとはPlamo Linuxでkernel 2.2ぐらいだったような. もう9年以上前か. 古いなあ. kernel 2.4→2.6の変更が一番大変だったかな(もちろんパッケージ管理なんてものは無い).
10月19日_
アルゴリズムはもういじれそうにないのであとは禁断のinline化しかないよなーとゴリゴリasm化してた. これ以上すると流石にコードキャッシュからはみ出るだろうというぐらいやってみると本当に遅くなった. 作為的な例じゃなくてリアル(?)なパターンで本当にそういうのに遭遇すると面白い. あとは適当にinlineしすぎた部分をcallに置き換えていく作業が必要だ.
10月15日_

初めてのSymPy

昨日と同じことを今度はSymPyでもやってみる. 一部うまくいかない(分からない)ところがあったけど一応答えは出た. SymPy歴30分なのでよりよい方法があれば教えてください.

SymPyはPythonのライブラリなのでPythonをインストールしてから入れる. pythonを起動して対話的に

from sympy import *

まず使いたい変数を指定する. Symbolかvarを使う. varだとまとめて定義できて便利そう.

x = Symbol('x')
t = Symbol('t')
a = Symbol('a')
b = Symbol('b')
c = Symbol('c')
d = Symbol('d')
e = Symbol('e')
# または
var('x t a b c d e')

Pythonでは巾乗は^ではなく**なのに気を付ける.

f=integrate((exp(x)-(a+b*x+c*x**2+d*x**3+e*x**4))**2,(x,0,t))

solveに渡す式は等式ではなく関数なのに注意. 関数の零点を求めるということ. 最初ちょっと悩んだ.

ans=solve([diff(f,a),diff(f,b),diff(f,c),diff(f,d),diff(f,e)],[a,b,c,d,e])

答えは辞書で返ってきてるのでans[a]とかで参照する. t = 0での極限を出してみる.

>>> limit(ans[a], t, 0)
1
>>> limit(ans[b], t, 0)
1
>>> limit(ans[c], t, 0)
1/2

うまく行ってる感じ. ところがans[d]やans[e]はエラーになった.

>>> limit(ans[d], t, 0)
  File "<stdin>", line 1, in <module>
  File "C:\p\Python\lib\site-packages\sympy\series\limits.py", line 116, in limit
  (...)
raise ValueError('(%s).series(%s, n=8) gave no terms.' % (f, x))
ValueError: ((-4*_w**4*exp(_w) + 5*_w**4 + 84*_w**3*exp(_w) + 96*_w**3 - 780*_w*
*2*exp(_w) + 840*_w**2 + 3720*_w*exp(_w) + 3840*_w - 7560*exp(_w) + 7560)/_w**8)
.series(_w, n=8) gave no terms.

テイラー展開か何かが最大n = 8で決め打ちなのかな. 何かで変更できるのかちょっと調べたけど分からなかったのでとりあえず進める. 近似値をえるにはevalf()を使う. が, ans[a]に対してどうevalすればよいのか分からなかったので代わりにlimitを使った. 分母が0にならないところの値なので大丈夫だろう.

>>> limit(ans[a],t,log(2)/4096).evalf(50)
1.0000000000000000000000045896459498312533391085579
>>> limit(ans[e],t,log(2)/4096).evalf(50)
0.041670192358256038810322743017020389259233306763901

さて, もとに戻って上記limitのエラーの部分を見てみる. 実装コードをすぐ見られるのは利点やね. sympy/series/gruntz.pyの508行目.

def calculate_series(e, x, skip_abs=False, logx=None):
    """ Calculates at least one term of the series of "e" in "x".

    This is a place that fails most often, so it is in its own function.
    """

    f = e
    for n in [1, 2, 4, 6, 8]:
        series = f.nseries(x, n=n, logx=logx)
        if not series.has(O):
            # The series expansion is locally exact.
            return series

        series = series.removeO()
        if series:
            if (not skip_abs) or series.has(x):
                break
    else:
        raise ValueError('(%s).series(%s, n=8) gave no terms.' % (f, x))
    return series

うむ, どうみても8までしかやってない. しかも数値固定. 試しにforループのところを10まで増やしてみる. お, 今度はうまく行った. なんでパラメータで変更できるようになってないのだろう.

>>> limit(ans[d],t,0)
1/6
>>> limit(ans[e],t,0)
1/24

全般的に速度はMaximaよりもっさりという感じ. Mapleは流石にそれらよりずっと速い(10倍は違うと思う). でもPythonなのは何かと使いやすいなあ. 最後に全体のコード

from sympy import *
var('x t a b c d e')
f=integrate((exp(x)-(a+b*x+c*x**2+d*x**3+e*x**4))**2,(x,0,t))
ans=solve([diff(f,a),diff(f,b),diff(f,c),diff(f,d),diff(f,e)],[a,b,c,d,e])
limit(ans[e],t,log(2)/4096).evalf(50)
limit(ans[e],t,0) ## コードをいじらないとエラー
10月14日_

初めてのwxMaxima

こないだのx86/x64最適化勉強会の私の発表[高速な倍精度指数関数expの実装]のおまけのところでやってた数式処理ってどうするのという質問があったので書いてみる. 私はMapleを使ってたのだけどMaximaはフリーソフトで同じことができる. なお私はMaxima歴30分なのであしからず.

今回やりたいことはt = log(2) / 4096が与えられたときに

f:=∫0t(exp(x) - (a + bx + cx2 + dx3 + ex4))2dx

を最小とするようなa, b, c, d, eは何か?という問題(expの計算ではb=1固定でやった)を手で計算せずにコンピュータにさせること. まずは積分させてみよう.

integrate((exp(x)-(a+b*x+c*x^2+d*x^3+e*x^4))^2,x,0,t);

と入力してみる. 関数をxについて区間[0, t]で積分せよということ. おっと入力の終わりは[Enter]ではなくctrl-[Enter]なので注意. すると

Is t positive, negative, or zero?

と聞いてくる. へぇ. positiveと打つと無事計算結果が出た. 最初からtが正を仮定しておくには

assume(t > 0);

としておけばよいようだ. 関数の定義は := を使う.

f(a,b,c,d,e):=integrate((exp(x)-(a+b*x+c*x^2+d*x^3+e*x^4))^2,x,0,t);

fの最小値を求めるにはfをa, b, c, d, eで微分してそれらを0と置いた連立方程式を解けばよい. 微分はdiffを使う. たとえばaで微分するには

diff(f(a,b,c,d,e),a);

でOK. 連立方程式を解くにはsolveを使う. 複数の式は[]でくくる. たとえば

solve([2*x+3*y=13,4*x+y=11],[x,y]);

とするとx = 2, y = 3と表示される. 今回欲しい値はtの関数なのでそれをans(t)と定義してみよう.

ans(t):=solve([diff(f(a,b,c,d,e),a)=0,
               diff(f(a,b,c,d,e),b)=0,
               diff(f(a,b,c,d,e),c)=0,
               diff(f(a,b,c,d,e),d)=0,
               diff(f(a,b,c,d,e),e)=0],[a,b,c,d,e]);

えらく長い式が出力された. 正しく計算できてるかt → 0での値を見てみよう. 極限値を取るにはlimitを使う.

limit(ans(t),t,0);
[[a=1,b=1,c=1/2,d=1/6,e=1/24]]

ちゃんとexp(x)をテイラー展開したときの最初の4項の係数がでている. よさそうだ.

数値の定義は:を使う.

t : log(2)/4096

こうしてからans(t)とすると厳密解が表示される. 今回欲しいのは小数の値なのでfloatを使ってみる.

float(ans(t));
[[a=-7482132.582755132,
  b=2.646673732478273*10^11,
  c=-7.036349605814977*10^15,
  d=6.479223200922161*10^19,
  e=-1.912918493656289*10^23]]

あれ, 値があきらかにおかしい. 実はMapleでも同じことが起こって少し悩んでいた. 精度を大きく取らないと誤差が大きくてよい近似値にならないようだ. 内部演算の精度を変更するにはfpprecに値を設定すればよい. 任意精度の計算にはbfloatを使う.

fpprec : 50;bfloat(ans(t));
[[a=1.0000000000000000000000045896015893185326124785262b0,
  b=9.9999999999999999918636653907610524181757630578724b-1,
  c=5.0000000000003365597465341327601927431869990904259b-1,
  d=1.6666666613631816947174247945213363211134669284956b-1,
  e=4.1670192348702792353452213368172890252007837484349b-2]]

今度はちゃんと出た. 今回の操作を全て書いてみよう.

assume(t>0);
f(a,b,c,d,e) := integrate((exp(x) - (a + b*x + c*x^2 + d*x^3 + e*x^4))^2, x, 0, t);
ans(t):=solve([diff(f(a,b,c,d,e),a)=0,
               diff(f(a,b,c,d,e),b)=0,
               diff(f(a,b,c,d,e),c)=0,
               diff(f(a,b,c,d,e),d)=0,
               diff(f(a,b,c,d,e),e)=0],[a,b,c,d,e]);
fpprec : 50;bfloat(ans(log(2)/4096));

これだけでよい. すばらしい. 今検索するとはじめてのMaximaという膨大な(1000ページ超!)Maximaと数学の解説書を書いてる方がいらっしゃった. ジャンルが多岐に渡っていて目次を見ているだけでため息が出る.

10月13日_
ICFP Programming Contest 2011 舞台裏&表
JRが止まってていくのに難儀した. ICFPのコンテストをどうやって作るのか, どうやって解くのかという話. ゲームバランスを決めるのが大変らしい. 今回は非常にうまくいったのでぜひやってみてとのこと. 懇親会にも参加したかったが, まだ体調が悪くて諦める.
・1割速くなったと思ったのは勘違い. まだあと1割速くせねば.
10月12日_
カーネル読書会のこと, 勉強会の問題の典型例
twitterでもいろいろみんなコメントしてるけど, 答えは明白とちゃうのか. hyoshiokさんが最初に
わたしのモチベーションが下がったというのもあるのだが(まあ、それが唯一の原因であるのだけど)
って書いてるやん. 開催メリットが低くなっちゃったらどうしようもない. いくら開催コストを下げるアイデアを出したとしても, さまざまな理由で却下されるだろう. ここは初心に戻って最新kernelのソースのどこかを読んでモチベーションを回復させるか, あるいは回復するまでしばらく放置でええんじゃないかな.
・他の人が書いたasmのコードの中に無駄なmovを発見. これで約1割速くなった. あとちょっと. しかし, 自分のじゃないasmコードのチェックはかなり大変だなあ. あんましやりたくない.
10月9日_
・アルゴリズム的にはぎりぎりな感じ. とするとレジスタ退避と復元回数を減らすしかないか.
・せっかくの休みなのに喉がひどく痛くて鼻水もでる. 早くちゃっちゃと終わらせたいのだけどなあ.
10月8日_
うーむ, まだ1割速くしないといけないのか. かなりきついな.
10月7日_
他人が書いたコードから不要なmemcpyを探してちまちまと削る. 慣れてないと意外とcriticalなところで入り込んでしまうんだよな.
10月1日_
x86/x64最適化勉強会2終了
みなさまお疲れさまでした. kzmogiさん, lunatic_starさん, (株)AXELL社の方々, 会場提供とustreamの配信どうもありがとうございました. とても快適な会議室でした.
とりあえず資料へのリンクとtogetter.
  1. maropu : 研究動向(Intel Lab.等)からみるx86/x64最適化手法[ust-1]
  2. kzmogi : 適応フィルタのSIMD最適化, x86opti_02_kzmogi.zip [ust-2]
  3. w_o : perfの使い方 [ustなし]
  4. herumi : JIT profiling with VTune [ustなし]
  5. t_teruya Kronecker-Jacobi記号の実装, github.com/tell/num : [ust-3]
  6. shinh : edit distanceの高速化 [ust-4]
  7. herumi : 高速な倍精度指数関数expの実装, fmath [ust-5]
  8. msiro : Intel Architecture Code Analyzerについて [ust-6]
次回はあるとしたら来年初頭かなあ. 発表したい方はぜひ連絡を!

御意見は光成滋生<herumi@nifty.com>までお願いします