JavaScript:float32による小数演算最適化

Posted on

また翻訳です!


コンピュータ上で浮動小数点を表現する方法はいくつか存在する。だいたいのアーキテクチャではIEEE754標準を使っている。この標準では、倍精度浮動小数点(a.k.a. double or float64)を64bitで表し、単精度浮動小数点(a.k.a. float or float32)を32bitで表している。名前が示す通り、float64はfloat32よりもっと正確なんだ。だからパフォーマンスがあんまり重要視されないときは大体こっちを使おうねってアドバイスされる。

でもね、float32を使うことにもこんな利点がある:

  • float32は計算操作をするのにfloat64よりCPUサイクルがしばしば少ない。精度が低いからだ。
  • float64の計算をfloat32に対して行うには、いくらかオーバーヘッドが存在する。なぜなら、一旦float64に変換しないといけないから。これは、Float32arrayを使ってデータを格納しているときにしばしば発生する。WebGLを使ったコードや、EmscriptenでコンパイルされてC++のfloatがfloat32として格納されいるけど計算はJavaScriptのfloat64として行っている場合ではよくあることだ。
  • float32は半分のメモリしか使わない。プログラムのメモリ総使用量を減らすし、メモリーがボトルネックになるような(memory-bound)計算では速度も速くなるだろう。
  • float32に特化した標準C関数は、我々の調査したいくつかのライブラリではfloat64のバージョンより速かった。

JavaScriptにおいて、Number値はfloat64として定義されていて、すべての計算もfloat64として行うことが定義されている。さて、JS処理系がいままで長い間高速化のために利用していたfloat64の特性が、float64の値が十分小さい整数のとき、float64での計算結果は整数での計算結果と全くおなじ、ということだ(桁溢れが発生しなければ!)。JS処理系は、この特性を、整数の値を(floatでなく)生の整数として保持しておいて、整数演算で計算することで(そして桁溢れをチェックすることで)利用してきた。実際、次の少なくとも3種類の表現が実際に使われている:31bit整数、int32_t、uint32_t(くわしくはAndy Wingoの書いたこの値の表現に関する記事を見てみてね!)。

これらを全部踏まえると、一ついい感じの疑問が湧いてくる:似たような事を、float32で出来ないかな!?

答えは「しばしば可能」で、ES6で追加される新しいMath.fround関数を使うとプログラマ側がいつfloat32演算を使うかをコントロールすることが出来る。

いつ安全にfloat64演算の代わりにfloat32演算が使える?

Samuel A. Figueroaが“When is double rounding innocuous?” (SIGNUM Newsl. 30, 3, July 1995, 21-26)で書いているんだけど、float32とfloat64で交換がなりたつ恒等式がいくつかある。入力がfloat32であるかぎり、float32での演算とfloat64での演算、どちらを使ってもよくて、同じ結果を得られるんだ。もっと正確に言うなら、{+,-,*,/}のそれぞれの演算について、float32版の演算を行う関数であるop_f32と、float64版の演算を行う関数であるop_f64があったとすると、次の恒等式が成り立つんだ。C++で表現すると、こんな感じ:

  float op_f32(float x, float y);
  double op_f32(double x, double y);
  assert(op_f32(x,y) == (float) op_f64( (double)x, (double)y ));

似たような恒等式が、sqrtとかその他の数学関数のような単項演算にも成り立つ。

ただし、この性質は、それぞれの演算の前後でfloat32にcastした場合に成り立っている。たとえば、x=1024でy=0.0001で、さらにz=1024であるときに、(x+y)+zはfloat32の二回の足し算として見た時とfloat64の二回の足し算とみた場合では、もはや同じ結果にはならない。

この恒等式は、コンパイラがfloat64の演算のかわりにfloat32の演算を使って最適化を行うための、前提条件を与える。実際に、gccはこの恒等式を利用していて、さっきの恒等式の右辺のような式(float64のコード)があったら、左辺のような式(float32の式)を生成する。

でも、いつJavaScriptはfloat32の値を持ちうるんだろう?WebGLとそれに付随するTypedArray APIが使えるHTML5では、答えはこうなる:Float32Arrayから読む時と、書く時、だ。

たとえば、このコードを見てみよう:

  var f32 = new Float32Array(1000);
  for(var i = 0; i < 1000; ++i)
    f32[i] = f32[i] + 1;

ループの中の足し算のコードは、たしかにさっきの恒等式と一致してる:f32配列からfloat32の値を持ってきて、float641に変換され、float64同士の足し算演算が行われ、そしてf32に書き戻すためにfloat32に再度キャストされる。(足している”1″というリテラルも、float32からfloat64にキャストしていると考えた。なぜなら、float32で1は正確に表現できるからだ)

でも、もっと複雑な式を書きたいときはどうしよう?不要なFloat32Arrayへの出し入れを式のそれぞれの式の一部分に追加して、計算した式の一部分の結果がつねにfloat32にキャストされるようにすることも出来る。でも、これらの追加の(不要な)出し入れはプログラムを遅くするし、我々の目標はプログラムを速くすることなので、それでは意味がない。

うん、十分に頭のいいコンパイラーさんなら確かに不要なこの手の出し入れを除去することができる。でも、性能の予測可能性は重要だから、もっとこの最適化が行われることを確実にしたい。かわりに、わたしたちは小さな組み込み関数を提案し、Ecma Script 6規格で取り入れられた:Math.froundだ。

Math.fround

Math.froundは来るべきES6標準で提案されている、新しい数学関数だ。この関数は入力された値をいちばん近いfloat32の値にまるめて、そのfloat32をNumber値として返す。つまり、Math.froundは意味としては2つぎのplyfill3と同値だ:

  if (!Math.fround) {
    Math.fround = (function() {
      var temp = new Float32Array(1);
      return function fround(x) {
        temp[0] = +x;
        return temp[0];
      }
    })();
  }

いくつかのブラウザではTypedArraysに対応していないけど、そういう時はまた別の複雑なpolyfillがあるよ。でも、TypedArraysを使わなくても、SpiderMonkey(FirefoxのJSエンジン)とJavaScriptCore(SafariのJSエンジン)の両方ですでにMath.froundが実装されているよ。V8(ChromeのJSエンジン)でも実装する計画があって、このissueでその状況が確認できる(訳注 すでに実装されてるっぽい)。

結果として、float32の演算として演算を行い続けるには、単に途中の結果をMath.froundで包めばいい:

  var f32 = new Float32Array(1000);
  for(var i = 0; i < 999; ++i)
    f32[i] = Math.fround(f32[i] + f32[i+1]) + 1;

プログラマに速いコードを書くのを許すだけではなく、EmscriptenみたいなJSコンパイラにも、float32のコードをよりbetterにすることを許す。たとえば、Emscriptenは現在C++のfloat演算をJSのNumberを用いたコードにコンパイルする。技術的には、Emscriptenは原理的にはFloat32Arrayへの出し入れを全ての演算ごとに行うことで余計なfloat64精度を取り除くことは出来るけど、とても遅くなってしまうだろう。だから、忠実性はパフォーマンスのために犠牲にされた。

この違いが何かを壊すことは非常にレアなんだけど(もし壊れてしまうなら、そのプログラムは「未定義動作」に依存してるということだ)、この違いがバグを起こした所を実際に見たことがあって、しかもそういうバグを追跡するのは全く楽しくない。Math.froundを使えば、Emscriptenはもっと効率的で、もっと元のC++コードに忠実なコードを生成できるんだ!

Float32 in IonMonkey

わたしのインターンシップでのプロジェクトは、Firefoxにこれらの最適化を組み込むことだった。まず最初のステップとして、一般的なFloat32の演算のサポートをIonMonkey JITバックエンドに追加した。つぎに、Math.froundを一般的な(最適化されていない)ビルトインとしてJavaScriptエンジンに追加した。最後に、Float32Array/Math.froundを使っているのを認識してfloat32演算を出来る限り使うように最適化する最適化パスをコンパイラに追加した。これらの最適化は、Firefox27から使うことが出来る(今はAuroraリリースチャンネルに入ってるよ→訳注 もうStableでもVersion30ですね)

で、どれくらい早くなるんだろう?ちいさなベンチマークでは(C++とJavaScript両方で)、50%くらいの大きなスピードアップを示した。でもこういうちいさなベンチマークはしばしばミスリーディングを引き起こすから、わたしは次のようなもっと現実的なベンチマークを書いて、float32の大量の演算がどれくらい速くなるか、実際にありそうな問題で調べた:

  • 逆行列計算:このベンチマークでは、いくつか行列を作って、その逆行列を計算し、float64とfloat32でどれくらい精度が失われるかを比べるために、もとの行列と掛け算した。このベンチマークには、gl-matrixの改造したバージョンを利用した。gl-matrixは、WebGLの行列計算を使ったフレームワークで、実世界のアプリケーションでも使われている。float32バージョンを作るためには、Math.froundの呼び出しだけを追加した。
  • 画像処理行列計算:このベンチマークでも行列をいくつか生成して、画像処理でよく使われる行列の計算を行う。移動・回転・スケーリング、などなど。このベンチではたくさんの基本操作や複雑な操作(例えばMath.cosやMath.sinを回転操作で呼び出す)をおこなう。で、 float32の演算を最適化したとき、このベンチではすごい改善が見られた。このベンチマークも、gl-matrixの改造したバージョンを使っている。
  • 指数関数:このベンチマークでは、大きなFloat32Arrayを予測可能な値で埋めて、つぎにそれぞれの要素の指数関数を、テイラー展開をつかって求めた。このベンチマークの主な目的は、足し算と掛け算の重さを測ることだ。
  • 高速離散フーリエ変換:このベンチマークでは、適当なニセのバッファを作って、高速離散フーリエ変換(FFT)を行った。このベンチマークには、他にも基本操作と数回のMath.sqrtへの呼び出しが含まれている。FFTのコードは、既存のライブラリ、dsp.jsから取ってきた。

次の表が違うデバイスで最新のFirefox Nightlyを実行したときの結果を示している。それぞれの数値が、Math.froundの呼び出しを追加してfloat32最適化を許すことで得られたスピードアップを示している(だから、高大きい方がより良いことを示している)。

使っているデスクトップマシンは、「ThinkPad Lenovo W530 (Intel(R) Core(TM) i7-3820QM CPU @ 2.70GHz, 8 cores, 16 GB RAM)」だ。電話もしくはタブレットと書いてあるときは、Firefox NightlyのAndroidバージョンの事をしめしている。

この結果をひと目みたら、自分でもベンチマークを取ってみたいって思うんじゃない?ここで自分で実行できるよ!(あ、でも忘れないで。Firefox 27かそれ以降を使うんだ)

ベンチマークのソースコードはgithubで見れるよ(gh-pagesブランチにある)

Device Matrix Inversions Matrix Graphics Exponential FFT
Desktop (x86) 33% 60% 30% 16%
Google Nexus 10 (ARM) 12% 38% 33% 25%
Google Nexus 4 (ARM) 42% 26% 38% 5%
Samsung Galaxy S3 (ARM) 38% 38% 24% 33%

Polyfilling Math.fround

すべてのJSエンジンでMath.froundが使えて最適化されるようになるまで、どうすればいいんだろう?上で紹介した忠実な方法ではなくて、こんな単純な恒等関数を使えばいい:

  var fround = Math.fround || function(x) { return x }

これがさっきのベンチマークが使っていて、かつ、上で示したように、殆どのコードで((訳註 一度他の関数に代入することによる?))違いをうまない。

すべてのモダンなJS処理系はだいたい小さな関数をハイエンドなJITでインライン化するので、このplyfillはパフォーマンス上のペナルティになることはない。それを確かめるために、例えばChromeの開発版で実行することもできる。しかしながら、いくつかの巨大なソースコードではインライン化が実行されず、plyfillのせいでパフォーマンスが悪くなったこともあった(たぶん、インライン化の深さ制限か、関数が何かしらの理由でJITコンパイルされなかったんだろう)。

まぁ簡単に言うと、「試してみる価値はある。そしてテストしよう」。

結論

結果は私達を元気づけてくれるものだった。Math.froundはES6の標準にはいるから、ほかのJS処理系も同じ最適化を行うようになってくれることが期待される。Webをゲームプラットフォームにしたいという観点からすると、低レイヤの今回のような最適化はより重要性を増してきていて、Webをよりネイティブ性能に近づけてくれると言えるだろう。ぜひFirefoxのNightlyAurolaでこの最適化をテストして、何かバグを見つけたらぜひ教えてね4

今回のことに参加してくれたすべての人に感謝したい:Jon CoppeardとDouglas CrosherにはARMバージョンを実装してくれたことを、Luke Wagner、 Alon ZakaiとDave Herman には校正とフィードバックを、そして、JavaScriptチームのみんなにヘルプとサポートを。


 

感想

JSの最適化もここまで来たかーって感じがします。今後はこういう「プログラマと処理系が手を組んで行う最適化」が増えてくるのでしょうかね。そういう最適化はたぶんまだ沢山できて効果がありそうな一方、プログラムを書いている人が知ってるか知ってないかでかなり速度が変わってくるわけで、JSを書くのがどんどん大変になっていきそうな…。やっぱJSはアセンブリなんや…(何)。

  1. 訳注 JavaScriptのNumberはfloat64と定義されているのでした []
  2. 訳註 速度を除いて []
  3. なんとWeb時代に生まれた単語で、古いブラウザで新しいブラウザと同じ機能を提供するための方法のこと []
  4. 訳註 もう30なので、製品版コードに入っています []

ドワンゴのC++勉強会でLTしてきたっ!

Posted on

土曜日に、ドワンゴC++勉強会 #1に行って発表してきましたっ!

タイムシフトがまだ見れます

最初、「コンパイルの話ばっかりなのかなぁ、わたし、そこまでコンパイル時処理は書いてないんだけど…」と思って普通に聞く側で行こうかと思っていました。

行きたいけど、でも定員の倍くらい既に埋まってて絶望的だなぁ…と思っていたら、「Cocos2dxの闇 @ponkotuy」というのが見えて、「そうかー、実行時の話してもいいんだ…!そういえば去年お仕事で無限にcocos2d-xとバトルしてpull requestとか送ってたなあ結構いろいろ…」と思ったので、その時に活躍したvalgrind(ヴァルちゃん)の話をしてLT枠で参加することにしました!ぶっちゃけcocos2d-xがC++製だったのを忘れていたのは秘密だ

valgrindも色々制限があったり(たとえばリークの検出は、ちゃんと最後まで実行できないとうまくいかない)、遅かったりと万能ではないのですが、それでも最終手段として使ってみるのはアリだと思います。その後に知ったのですがLLVMにも同じようにメモリのアクセスをチェックしてくれる機能があるらしいです。

cocos2d-xは…もう何も言うまい。UnityもUnityで☆闇☆ですけどね…。


格好は五月祭のときにつくったコンちゃんの衣装で してみました。コンパイル時の話が多かったので(?)

家からドワンゴに行くまでもこの格好で行ったのですが、ガン見する人とさらっと流す人が明確に分かれてたので面白かったです。ちなみに人からジロジロ見られる率ですが、

白衣>>>(職質の差)>>>コスプレ>>>>>>>>>>>>> 異性装

みたいな感じです。白衣着てるのがこの社会では一番変な人扱い、らしい。


その後は妖怪ハウスで江添さんのピザを頂いたりしました。シーフードピザ、美味しかったです!ありがとうございました!(小学生の日記みたいな文章だ)


Fortranが書きたくないなら機械語を埋め込めばいいじゃないっ!

Posted on

Fortranが書きたくない!!

わたしは大学で天気とか地震とか海流といった「地球物理」と呼ばれるジャンルを学んでいるのですが、このジャンルではコンピュータをガンガン回して計算しまくります。気象庁が毎日やってる天気予報もそうですし、「地球シミュレータ」でエルニーニョや地球温暖化の予測とかをやっているのもそうです。

で、そういう学科なので、コンピュータの実習があるのですが、なぜか今でもFortranを教えています。

たしかにFortranには、LAPACKのような優れた数値計算ライブラリがあったり、スパコンで書く時はコンパイラのサポートが充実してたりします。なので、そういう特定の用途には悪くない言語だと思うのですが、デスクトップPCで十分扱えるようなデータサイズで、そもそも計算速度はさして重要ではなくてデータ解析手法を学ぶ実習で使うのには辛すぎます。陽には書きたくありません。1でも書かないと単位くれないんだってさ!

なので、Fortranを書かずにFortranを書く方法を検討して、実際に実装しました。今回のソースコードも、すべてgithub上に上げています。

Fortranを書かずにFortranを書くには?

最初、Fortran上で小さなVMのインタプリタを実装して、そのVMの仮想機械語にコンパイルされる言語のコンパイラを実装しようと思いました。

が、VMを書くのが普通にしんどいのと、言語とコンパイラをこのためだけに実装するのもかなりつらいので、もっとラクな方法は無いか考えてみました。

冷静に考えると、別にVMのインタプリタを実装しなくとも、コンパイルされたFortranのコードは既にCPUという名前のx86機械語のインタプリタ機械(?)で実行されています。このインタプリタ(?)を利用できないでしょうか。利用できれば自力でVMのインタプリタを実装しなくていいので手間を省略できますし、コンパイラもgccとかclang/llvmとかすごいやつがいっぱい使えます!

具体的にはそのためにどうすれば良いのかというと、Fortranのinteger配列にC言語のソースをコンパイルしたx86の機械語の数字を沢山並べて入れて、そのinteger配列のポインタをC言語の関数ポインタに無理やりキャストして、そして実行すればいいのですっ!!

module a;interface;subroutine f (z) bind(c);use iso_c_binding;type(c_ptr), value :: z;
end subroutine;end interface;type G;procedure(f), pointer, nopass :: x;end type;type J;integer, pointer :: x;
end type;end module;program b;use iso_c_binding;use a;implicit none;type (J), pointer :: x(:);
real(8),pointer :: dbl(:);type(G), pointer :: p;integer, pointer :: d(:);integer :: i;integer :: n;allocate( x(1) );
allocate( dbl(10000000) );allocate( d(10000000) );dbl = 0;x(1)%x => d(645);d(1)=-443987883;d(2)=1213580125;d(3)=267576713;
d(4)=1223181585;d(5)=1223181709;d(6)=48087689;d(7)=450755289;d(8)=-128612024;d(9)=-398095544;d(10)=-532313784;
d(11)=1158680562;d(12)=1438866912;d(13)=1223002440;d(14)=1223705997;d(15)=-338050423;d(16)=-1991763235;
d(17)=-1958152107;d(18)=-1991708603;d(19)=267577413;d(20)=1575503120;d(21)=-1991748157;d(22)=286257893;
d(23)=-1924601787;d(24)=-1991710651;d(25)=-654123582;d(26)=1209720318;d(27)=1224234377;d(28)=1223181707;
d(29)=-220183159;d(30)=-532344817;d(31)=1213580125;d(32)=267576713;d(33)=1223181585;d(34)=1223181709;d(35)=48087689;
d(36)=450756569;d(37)=-128612024;d(38)=-398095544;d(39)=-532313784;d(40)=1158680562;d(41)=1438866912;
d(42)=-219838136;d(43)=-398126833;d(44)=-398095032;d(45)=-574453432;d(46)=-572401406;d(47)=1435060250;
d(48)=1166756088;d(49)=1166625000;d(50)=269480672;d(51)=-1017257915;d(52)=-443987883;d(53)=-394426040;
d(54)=-1192987255;d(55)=0;d(56)=-129660600;d(57)=16008647;d(58)=-352321536;d(59)=-196768982;d(60)=-1924622264;
d(61)=50452;d(62)=-1958215680;d(63)=21555269;d(64)=269480656;d(65)=269480448;d(66)=267581517;d(67)=267567448;
d(68)=-2080881391;d(69)=-1962806203;d(70)=1161557061;d(71)=1221491940;d(72)=1224230283;d(73)=-220707447;
d(74)=-666562545;d(75)=1213580125;d(76)=-1991711351;d(77)=-1991710595;d(78)=1435099253;d(79)=47324;d(80)=-1991770112;
d(81)=1170733125;d(82)=244;d(83)=-1958286592;d(84)=-1740049339;d(85)=-988508856;d(86)=0;d(87)=-398095544;
d(88)=-221249208;d(89)=-1962405873;d(90)=-1740049339;d(91)=-988508856;d(92)=0;d(93)=-532313272;d(94)=-221249208;
d(95)=-234876913;d(96)=-222209777;d(97)=-129167345;d(98)=-1051193358;d(99)=1158746098;d(100)=-196770824;
d(101)=-196769023;d(102)=2094810427;d(103)=1166756018;d(104)=1166625016;d(105)=269480656;d(106)=-1017262011;

(略)

d(761)=-800748728;d(762)=-389576376;d(763)=-1168;d(764)=1223710091;d(765)=1220572555;d(766)=-1178057333;d(767)=20;
d(768)=-389576376;d(769)=-683;d(770)=-1980742261;d(771)=535478722;d(772)=-120467455;d(773)=-223591031;
d(774)=-1404753393;d(775)=-2008708280;d(776)=1118194;d(777)=16008647;d(778)=-352321536;d(779)=-196768970;
d(780)=-2092394424;d(781)=-1924660800;d(782)=50452;d(783)=-1958215680;d(784)=21530693;d(785)=-196768830;
d(786)=-1924622264;d(787)=50444;d(788)=-1958215680;d(789)=21545029;d(790)=9128136;d(791)=-2096985784;
d(792)=-1962806203;d(793)=1161557061;d(794)=-1195213652;d(795)=0;d(796)=1213580233;d(797)=-1017256567;
call c_f_pointer(C_LOC( x(1) ), p);i=0;read (*,*), n;dbl(1) = real(n);do i=2,n+1;read (*,*), dbl(i);end do;call p%x( c_loc(dbl(1)) );
n = int(dbl(1));do i=2, n+1;write (*,*), i, dbl(i);end do;deallocate ( d );deallocate ( x );deallocate ( dbl );
contains;end program;

Fortran Side

Fortranでのキャストテクニック

今回はintegerの配列のポインタをCの関数へのポインタだと思わせたいため、(Cで言うと)int*をvoid (*ptr)(double*)のポインタにキャストする必要があります。しかし、Fortranは割りと型に厳しく、Cのように簡単にキャストしたり出来ません。なので、ちょっと撚(ひね)る必要があります。

!! 関数の型定義。
!! typedef void(*f)(void*);
interface
	subroutine f (z) bind(c)
		use iso_c_binding
		type(c_ptr), value :: z
	end subroutine
end interface
!! struct FPointer { f fptr; };
type FPointer
	procedure(f), pointer, nopass :: fptr
end type
!! struct IPointer { int* iptr; };
type IPointer
	integer, pointer :: iptr
end type

このように2つ構造体(IPointer,FPointer)を作っておき、この2つの構造体のポインタをC言語のポインタを経由することで無理やりキャストします。

type (IPointer), pointer :: ip(:);       !! IPointer* ip;
real(8),pointer :: dbl(:);               !! double* dbl;
type(FPointer), pointer :: fp            !! FPointer* fp;
type(C_PTR) :: cptr                      !! void * cptr;
integer, pointer :: d(:) !命令列用の配列  !! int* d;
allocate( ip(1) )                        !! ip = malloc(sizeof(IPointer));
allocate( dbl(10000000) )                !! dbl = malloc(sizeof(double)*10000000);
allocate( d(10000000) )                  !! d = malloc(sizeof(int)*10000000);
!! 命令列をセット
d(1)=-443987883;
d(2)=1213580125;
!! 略
d(790)=-443987883;
d(791)=50013;
!! 645番目の要素がC言語側エントリポイント関数の頭なので、そこへのアドレスをセット。
ip(1)%iptr => d(645);
!! IPointerへのポインタを、C言語のポインタ、いわばvoid*に落とす
cptr = C_LOC( ip(1) );
!! そこからさらにFPointerへのポインタにキャストする(専用の関数を使います)
call c_f_pointer(cptr, fp)
!! 無理やり変換した関数ポインタ経由でバイナリコードにジャンプ
call fp%fptr( c_loc(dbl(1)) )

IPointer 構造体へのポインタを一度Cのポインタ(いわばvoid*)にキャストし、それを再度void*からFPointerのポインタにキャストしています。構造体の中身の型が「intへのポインタ」と「関数へのポインタ」で違うので、構造体を通して、その中身についてキャストすることができました。

なんでこんな七面倒臭いことをしているかというと、intのポインタと関数のポインタを直接キャストすることができないからです。intのポインタは type(C_PTR)、関数のポインタはtype(C_FUNPTR)なので互換性がなく、さらにC言語みたいにintに無理やり落としてもう一度ポイ ンタにするような操作もできません。しかし、一回構造体でくるんでしまえばどちらも値のポインタとして扱い、キャストすることができます!

さらに、現在のLinuxでは、通常データ領域をプログラムとして実行することはできない(DEP)ので、コンパイルする時にその制限をはずす必要もあります。これには、 -z execstackというオプションスイッチが使えます。

これで、任意のx86機械語の列をFortranから実行することがきました。あとは、その機械語の列をどうやって作るかという話になりますね。原理的には自分で書いてもよいのですが、一応課題をやっているコードなので、結構複雑です。なので、C言語で書いたコードをコンパイルして、その結果を利用したいところです。

Binary Side

コードのコンパイルにはGCCを使うことにしましょう。

この時問題になるのは、Fortranの動的に配置した配列の上に機械語を置くので、機械語の置かれるアドレスが実行時になるまでわからないことです。なので、機械語がどこに置かれるかがわかっていることが前提の通常のコードはうまく動いてくれません。

そのようなケースは実は珍しくなく、共有ライブラリもアドレス上のどこにロードされるかは分からないので、ちゃんとそのためのコードをコンパイラが出力します2。そのようなコードのことを、PIC(Position Indepentent Code)と呼びます!gccでは、-fPICというコンパイルオプションをつけると、PICとしてコンパイルしてくれます。

じゃ、-fPICを付けてコンパイルしたバイナリから機械語をもってくれば、目的は達成できるのでしょうか…?

Position Independent Code

ところがどっこい、これではうまく行きません。なんでかというと、一個前に翻訳した記事で解説されている、PLTとGOTがあるからです。

こんな非常に簡単なC言語コードを考えてみましょう。

void calledFunction(){
}

void function(){
  calledFunction();
}

これをコンパイルするためのMakefileはこんな感じで、

.PHONY: all

all:
  gcc -c -o test.o test.c -fPIC
  ld -shared -fPIC -o test.so test.o

コンパイルした結果が、こんな感じでーす。

test.so:     file format elf64-x86-64


Disassembly of section .plt:

0000000000000280 <calledFunction@plt-0x10>:
280:    ff 35 82 0d 20 00        push   QWORD PTR [rip+0x200d82]        # 201008 <_GLOBAL_OFFSET_TABLE_+0x8>
286:    ff 25 84 0d 20 00        jmp    QWORD PTR [rip+0x200d84]        # 201010 <_GLOBAL_OFFSET_TABLE_+0x10>
28c:    0f 1f 40 00              nop    DWORD PTR [rax+0x0]

0000000000000290 <calledFunction@plt>:
290:    ff 25 82 0d 20 00        jmp    QWORD PTR [rip+0x200d82]        # 201018 <_GLOBAL_OFFSET_TABLE_+0x18>
296:    68 00 00 00 00           push   0x0
29b:    e9 e0 ff ff ff           jmp    280 <calledFunction@plt-0x10>

Disassembly of section .text:

00000000000002a0 <calledFunction>:
2a0:    55                       push   rbp
2a1:    48 89 e5                 mov    rbp,rsp
2a4:    5d                       pop    rbp
2a5:    c3                       ret

00000000000002a6 <function>:
2a6:    55                       push   rbp
2a7:    48 89 e5                 mov    rbp,rsp
2aa:    b8 00 00 00 00           mov    eax,0x0
2af:    e8 dc ff ff ff           call   290 <calledFunction@plt>
2b4:    5d                       pop    rbp
2b5:    c3                       ret

こんな感じで、ばっちりpltを経由して実行してしまいます。この機械語をFortran上に載せて強制的に呼び出したとしても、PLTがちゃんと解決されないとうまく動いてくれません。デフォルトのPLT解決処理は共有ライブラリとして実行されていることが前提になっているので、今回は利用することはできません。

integer配列のアドレスから関数のアドレスを解決することは原理的には可能なので、Fortran上でPLTのアドレスの所を埋めてちゃんと実行できるようにすること自体は可能(なはず)です。

しかし、かなり意味もなく複雑になってしまうので、できればPLTにジャンプするのではなくて、呼びたい関数に直接ジャンプするようなコードを出力して欲しいです。PLTを経由して関数を呼び出す時に使ってる0xe8命令は相対アドレスによるcall命令なので、可能なはずです。

リンカースクリプトを自分で書けばPLTとGOTを使わない…!?

リンク時にPLTを経由するようにリンクしているので、じゃあすごく単純なリンカースクリプト自分で書けばPLTみたいな便利機能使わなくなるんじゃない?

ということで、実際に書いてみました。こうなりました。

OUTPUT_FORMAT("elf64-x86-64", "elf64-x86-64","elf64-x86-64")
OUTPUT_ARCH(i386:x86-64)

ENTRY(ep);
MEMORY
{
ROM(rxai) : ORIGIN = 0, LENGTH = 64k
}
SECTIONS
{
.text : {} > ROM
.rodata : {} > ROM
.data : {} > ROM
. = ALIGN(4);
__bss_start = .
; .bss : {} > ROM
__bss_end = . ;
}

Makefileは、

linkerscript:
    gcc -c -o test.o test.c -fPIC 
    ld test.so test.o --script ./linker.script

こんな感じにすると指定できます。結果は、

test-with-linkerscript:     file format elf64-x86-64


Disassembly of section .text:

0000000000000000 <calledFunction>:
   0:    55                       push   rbp
   1:    48 89 e5                 mov    rbp,rsp
   4:    5d                       pop    rbp
   5:    c3                       ret    

0000000000000006 <function>:
   6:    55                       push   rbp
   7:    48 89 e5                 mov    rbp,rsp
   a:    b8 00 00 00 00           mov    eax,0x0
   f:    e8 ec ff ff ff           call   0 <calledFunction>
  14:    5d                       pop    rbp
  15:    c3                       ret

となって、見事、PLTを使わずに直接相対ジャンプするようになりました!

実は-pieだけで大丈夫

実をいうと、リンカースクリプトなんて書かなくても、-pieというオプションを付けて、通常の実行ファイルとしてコンパイルすればplt/gotを使わないコードを生成してくれます。

pie:
    gcc -c -o test.o test.c -pie
    ld -pie test.so test.o

-pieというオプションは、gcc –helpすると、

  -pie                     Create a position independent executable

ということで、実行時にどのアドレスに配置されても実行できる実行ファイルを生成してくれます。

実行ファイルは共有オブジェクトとは違って外部の他のプログラムからその中の関数が呼ばれるようなことはないので、内部のそれぞれの関数のアドレスを動的に解決する必要はありません。なので、PLTを経由せずに直接callするようなコードにしてくれる…のでしょう。

…でも冷静に考えると、共有ライブラリの中でならPLT経由せずに直接ジャンプでよくない?

全てが終わったあとに気づいたので時既に遅しという感じでしたが、一応書いておきます。

プログラムセクションを抽出する

残るところは、出来上がったファイルから機械語を抽出して、数字の列にするだけです。そのために、機械語が書いてある部分はファイルのなかでどこなのかを調べましょう。

セクションの一覧を出力するには、objdump -hが使えます。

% objdump -h test.so

test.so:     file format elf64-x86-64

Sections:
Idx Name          Size      VMA               LMA               File off  Algn
0 .interp       0000000f  00000000000001c8  00000000000001c8  000001c8  2**0
CONTENTS, ALLOC, LOAD, READONLY, DATA
1 .hash         00000028  00000000000001d8  00000000000001d8  000001d8  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
2 .dynsym       00000078  0000000000000200  0000000000000200  00000200  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
3 .dynstr       00000019  0000000000000278  0000000000000278  00000278  2**0
CONTENTS, ALLOC, LOAD, READONLY, DATA
4 .text         00000016  0000000000000294  0000000000000294  00000294  2**2
CONTENTS, ALLOC, LOAD, READONLY, CODE
5 .eh_frame     00000058  00000000000002b0  00000000000002b0  000002b0  2**3
CONTENTS, ALLOC, LOAD, READONLY, DATA
6 .dynamic      000000c0  0000000000200f40  0000000000200f40  00000f40  2**3
CONTENTS, ALLOC, LOAD, DATA
7 .comment      0000002c  0000000000000000  0000000000000000  00001000  2**0
CONTENTS, READONLY

機械語本体が置かれているのは、.textというセクションです。この中のデータをFortranに書きましょう。File offの項目が、ファイルの先頭からの位置で、Sizeという項目がセクションのサイズで、今回の場合、要はプログラムコード全体のサイズです。

さらに、関数呼び出しするには、呼び出したい関数の、プログラムコード内での相対位置も必要です。これには、nmコマンドが使えます。

% nm test.so
0000000000201000 D __bss_start
0000000000000294 T calledFunction
0000000000200f40 d _DYNAMIC
0000000000201000 D _edata
0000000000201000 D _end
000000000000029a T function
0000000000000000 d _GLOBAL_OFFSET_TABLE_
                 U _start

セクションの所のVMAもしくはLMAというのがメモリ上での位置なので、機械語列の中での相対位置を調べるにはここでのアドレスから.textセクションのVMA/LMAの値を引けば求めることができます。

シェルスクリプトを使うともちろんこの作業は自動化できます。呼び出す関数名は”ep”としました。EntryPointの略ですね!

_start=$(objdump -h test.so | grep \\s.text\\s | awk '{ print $6 }')
_size=$(objdump -h test.so | grep \\s.text\\s | awk '{ print $3 }')
_ep=$(nm test.so | grep \\sep$ | awk '{ print $1 }')

Cのソースをコンパイルする所から、命令列を抽出してFortranのソースに変換する作業まで、これらの作業はすべて自動化しているので、よかったらgithubのソースを見てね!

今回書き込んだinteger型は4バイトなので、エントリポイントのオフセットが4の倍数でないとき(よくある)には、その整列条件を満たすために数バイトズラしたりしないといけません。そのコードとかも入ってます。

Programming Side

これで終わり…ではありません!もう一個、プログラムを書くという問題があります!!今回の環境では、

  • 外部にあるライブラリ関数などは呼び出すことができない

という、非常に大きな制限があります。実行時に決定されるライブラリ関数の実際のアドレスを知らないし、知るためにはライブラリ関数を呼び出さないといけないので、「服を買う服がない」状態なんです。もちろん、Fotranでなら調べられるので、調べた結果をC言語の関数に投げればいいんですけど、Fortran書きたくないし(ぉ)。

科学計算すればいいだけなので、sin/cosとsqrtくらいが使えればなんとかなります。sin/cos/sqrt/absは、FPU命令を使って計算できるので、ライブラリは必要ありません。

そのためには、インラインアセンブリで陽にFPU命令を実行する必要があります。

double fsin(double v)
{
	intptr_t d0;
	__asm__ (
		".intel_syntax noprefix;"
		"fld qword ptr [%0];"
		"fsin;"
		"fstp qword ptr [%0];"
		".att_syntax;"
		: "=&r"(d0)
		: "0"((intptr_t) &v)
		: "eax", "ecx"
	);
	return v;
}

こんな感じのを、FPUの命令ごとに実装していきます。

関数だけでなく、piもFPUを使って取得しています。なんでかというと、今回は機械語だけを埋め込んで、それとは別にある定数の領域をFortranのソースに含めていないからです!

Cのソースに円周率の値を直接書いていてもすこし複雑な処理になると定数領域に飛ばされてしまって、.textセクションに含まれなくなってしまい、結果としてプログラムがおかしくなってしまいました。定数領域もFortranに埋め込むだけでその問題は回避できるのですが、その改修がちょっとめんどくさかったのでやめて、FPU命令を毎回呼ぶことで回避しました。

入力と出力ですが、実際のデータの読み込みと結果の出力のみ、Fortranが担当することにしました。C言語は関数のシグネチャをvoid (*fun)(double *)とし、doubleの配列を処理してその結果をもとの配列の領域に書き込むことだけに専念していただきます。

他の言語でやってみても面白いかもね!

と思った。Fortranはなんだかんだ任意の機械語を実行する難易度は低い言語だったと思います。

ここまでやるならFortranで書いたほうが早くね?

ちなみに

書かなくて良くなった(╹◡╹)

ひさびさに色々バイナリが触れてたのしかったです。まぁ結果としては良かったのではないでしょうか!

  1. っていうか、いままでプログラム書いたことがない人にいきなりこれ教えるのってどうなんでしょう…。 []
  2. Windowsでは、64bit版からそうなりました。32bitでは固定だそうです。 []

ELFの再配置シンボルの解決

Posted on

この記事は、

Resolving ELF Relocation Name / Symbols

の翻訳です。認めてくれたLeaf SRの人ありがとう!訳はがんばりましたが、間違ってる所もあるかもしれません。そこはご了承ください…。


共有オブジェクトのELFファイル内関数へのcall命令(たとえば、puts関数の呼び出しとか)は、直接関数のアドレスへ飛ぶのではなくて、PLT(Procedure Linkage Table)に飛ぶよ。PLTを使うと、関数の実際のアドレスを実行時に解決することができる。言い換えると、共有オブジェクトが実際にどこにロードされるかは実行時にならないと分からないので、PLTを使って実際のアドレスを解決している。

次のELFのサンプルを見てみよう:

ここに、ELFのテキスト・セグメント1内にある、アドレス0x804833cへの呼び出し命令がある:

$ objdump -d example-elf | grep 804833c | grep call

804843f: e8 f8 fe ff ff call 804833c

それで、次がこのサンプルファイルのPLTだ。さっきの命令で呼んでたアドレス0x804833cに注意してね。ここには実際には単なる*0×8049684へのジャンプ命令が並んでいる。

$ objdump -d example-elf | grep “section .plt:” -A 31

Disassembly of section .plt:

080482fc <__gmon_start__@plt-0x10>:
80482fc:       ff 35 70 96 04 08       pushl  0x8049670
8048302:       ff 25 74 96 04 08       jmp    *0x8049674
8048308:       00 00                   add    %al,(%eax)
     ...

0804830c <__gmon_start__@plt>:
804830c:       ff 25 78 96 04 08       jmp    *0x8049678
8048312:       68 00 00 00 00          push   $0x0
8048317:       e9 e0 ff ff ff          jmp    80482fc <_init+0x18>

0804831c <__libc_start_main@plt>:
804831c:       ff 25 7c 96 04 08       jmp    *0x804967c
8048322:       68 08 00 00 00          push   $0x8
8048327:       e9 d0 ff ff ff          jmp    80482fc <_init+0x18>

0804832c <__stack_chk_fail@plt>:
804832c:       ff 25 80 96 04 08       jmp    *0x8049680
8048332:       68 10 00 00 00          push   $0x10
8048337:       e9 c0 ff ff ff          jmp    80482fc <_init+0x18>

0804833c :
804833c:       ff 25 84 96 04 08       jmp    *0x8049684
8048342:       68 18 00 00 00          push   $0x18
8048347:       e9 b0 ff ff ff          jmp    80482fc <_init+0x18>

0804834c :
804834c:       ff 25 88 96 04 08       jmp    *0x8049688
8048352:       68 20 00 00 00          push   $0x20
8048357:       e9 a0 ff ff ff          jmp    80482fc <_init+0x18>

(注意:*は、C言語のポインタ参照と同じで、0×8049684に書いてある値、の意味)

*0×8049684が実際には何なのかを知りたければ、Global Offset Table (GOT)を探せばいい。こんな感じだ:

$ objdump -s example-elf | grep got.plt -A3

Contents of section .got.plt:
804966c 98950408 00000000 00000000 12830408  ................
804967c 22830408 32830408 42830408 52830408  "...2...B...R..

0×8049684には、42830408(リトルエンディアンで0×08048342)が書いてあった。ここでPLTに戻ってもう一回0×08048342を見てみると、jmp 命令でこのアドレスに飛んだあと、”push $0×18″という命令を実行することが分かる。0×18=24っていうのは、再配置テーブルへのオフセットアドレスだ。このpush命令に続いて、PLTの最初のアドレスへのジャンプである、”jmp 80482fc”が実行される。

この80482fcに書いてある命令列は、PLTの他の部分と大分ちがっているのはすぐ分かるだろう。最初の2つの命令、”pushl 0×8049670″と”jmp *0×8049674″は結構重要だ。

最初push命令はアドレス0×8049670をスタックにpushしていて、これはGOTを指している。次の(*0×8049674)へのジャンプだけど、このアドレスもGOT内のアドレスだ。この2つのアドレスは、ELFファイルの中ではどちらも”0″が書いてある。

というのも、これらは実行時に動的に埋められるからだ。実行時には、最初のアドレスの値は特定のライブラリが使われているかを識別するための番号になり、次のアドレスにはリンカーのシンボル解決ルーチンのアドレスが入る。これらのルーチンは前にスタックにpushされた0×18を、正しい場所を解決するために使う。

このテクニックは「遅延リンク(lazy linking)」と呼ばれている。再配置された関数のアドレス解決が、実行時に、しかも必要な時に、実際に関数が呼ばれた時だけ行われるからだ。

一度リンカーによって解決されると、リンカーはGOTエントリを書き換えて、最初の jmp *(アドレス)が、それまでのpush $0×18をするコードのアドレスではなく、解決された関数アドレスに直接ジャンプするようになる。

それじゃあ、どうやって関数のアドレスを解決するためのシンボル名を取ってくるんだろう?関数のアドレスを解決するには、’0×8049688′ではなくて、’sprintf’という文字列が必要だ。

SHT_RELというタイプを持っているセクションには、次のような構造体が並んでいる:

 typedef struct
 {
      Elf32_Addr r_offset;    /* Address */
      Elf32_Word r_info;    /* Relocation type and symbol index */
 } Elf32_Rel;

さて、readelfでELFファイル内の再配置エントリを見てみよう2

$ readelf -r /testbins/sha1

 Relocation section '.rel.dyn' at offset 0x420 contains 3 entries:
 Offset     Info    Type            Sym.Value  Sym. Name
 0804b54c  00001106 R_386_GLOB_DAT    00000000   __gmon_start__
 0804b598  00000505 R_386_COPY        0804b598   stderr
 0804b59c  00000d05 R_386_COPY        0804b59c   stdin

 Relocation section '.rel.plt' at offset 0x438 contains 2 entries:
 Offset     Info    Type            Sym.Value  Sym. Name
 0804b55c  00000107 R_386_JUMP_SLOT   00000000   feof
 0804b560  00000207 R_386_JUMP_SLOT   00000000   putchar

これらの情報はいったい何処から来たんだろう!?Offset/Info/Type/Value/Nameのそれぞれの情報、さっきの構造体に無かったよね?えっとね、これらの値のうち、殆どは実はr_infoフィールドから直接もってこれるか、r_infoフィールドを間接的に使うともってこれるんだ。

offsetの値はPLTの中を指していて、ここに再配置エントリの実体が置かれる。バイナリを逆アセンブルすると、PLTの中の正しいオフセットにある、このアドレスへjmpする命令が見えるはずだ。

今注目したいのは、r_infoフィールドだ。このフィールドに関係する、2つの(ほんとの事をいうと3つあって、3つ目は1と2を組み合わせると得られる)重要なマクロがelf.hに存在する:

#define ELF32_R_SYM(val) ( (val) >> 8)
#define ELF32_R_TYPE(val) ( (val) & 0xff)

これらのマクロはr_infoからそれぞれ別の値を得るためのマクロだ。いっこずつ見てみよう。SHT_RELタイプをもつセクション(さっきからみてる、再配置セクション)は他にsh_linkっていうメンバを持っている。このsh_linkは重要で、これらの再配置情報のシンボル情報(関数名とかね)を持つ、ほかのセクションを指しているんだ。このsh_linkは、gccでは大抵「dynsym」というセクションを指している。

/bin/lsのdynsymを実際に呼んでみた結果がこれ(読みやすくするために、いくつか省略してるよ):

$ readelf -S /bin/ls

   There are 26 section headers, starting at offset 0x126d8:

 Section Headers:
 [Nr] Name Type Addr Off Size ES Flg Lk Inf Al
 ...
 [ 4] .dynsym DYNSYM 080484a0 0004a0 0006b0 10 A 5 1 4
 ...
 [ 8] .rel.dyn REL 08049170 001170 000028 08 A 4 0 4
 [ 9] .rel.plt REL 08049198 001198 0002f0 08 A 4 11 4
 ...

セクション7と8のsh_linkメンバを見て。それぞれの4が、dynsym(dynamic symbol table)セクションを指している。この第4セクションが、再配置エントリに対応するシンボルネームが書いてあるセクションだ。ここまでやってきたことをまとめておくと、(1)SHT_RELタイプがついてるセクションを探して、(2)sh_linkメンバーの指し示すセクションをたどった。

よし、色々数字があるのは分かったから、あとはそれをどう組み合わせていけばいいんだろう。再配置テーブルのエントリを一個一個見ていって、r_infoフィールドの値を持ってきて、ELF32_R_SYM(val)マクロを使って値を取り出していけばいい。この数字はdynsymセクション(か、sh_linkの指している他のセクション)内のエントリに対応している。で、dynsymテーブルをパースしてエントリを調べれば、シンボルネームが解決できる。

他のマクロは何なんだろう?ELF32_R_TYPE(val)は再配置がどんなタイプなのかを教えてくれる。その値はelf.hにR_386_GOT32とかR_386_JMP_SLOTみたいなのが定義されている。この定義を使うと、再配置テーブルが関数のためのエントリなのかどうかとかが分かる(でも、rel.pltセクションの中のはだいたいそうだ3)。

ただし、ちょっと注意。全てのELFファイルがセクションヘッダを持っているわけじゃない。そういう時は、プログラムヘッダを使うと、動的なセグメントがどこかわかるし、再配置テーブルのアドレスとか、再配置テーブルのシンボルテーブルがどこかとかが分かる。

もしこのエントリに間違いがあったらぜひ教えてね!

  1. 訳注:テキストセグメントにはプログラムのコード本体が置かれている []
  2. 訳註:これから読むELFは、さっきまでのとは違うELFファイルだから、PLTのアドレスとかは違う []
  3. 訳註:PLTはProcedure Linkage Tableの略なのです []

洋裁経験ゼロから2日でつくるコスプレ衣装

Posted on

いろいろあって大学5年目ですが、先週末の学園祭には元気に参加しました!

今年の学園祭では、好きすぎて京都まで聖地巡礼に行ってきたアニメ、「いなり、こんこん、恋いろは。」のコンちゃんのコスプレをしてみました!

SONY DSC

SONY DSCしっぽも生えてるんですが、この写真だと見えないですね~。この格好でお客さんにエルニーニョの説明をしました。あんまりうまく説明出来ませんでしたけど!

たぶん「いなこん」を全話見た皆様でも「こんなキャラ居たっけ…?」みたいな感じかと思われますが、OPで一瞬だけ出てきます。

20140520_2

いなりちゃんといつも一緒にいる子狐の神使・コンちゃんが人間の姿で顕現する時の姿で、漫画ではみんなで一緒に海に行く話のあとの回で登場します。

20140520_3

何度見てもコンちゃんかわいい~!!(フサ フサ)

この話は殆ど他の話のフラグに依存しないし影響もしない、いわゆる「エピソード」なので、10話しかないアニメではカットされたのだろうと思います。漫画8巻の付属BDにはこのエピソードのアニメが入るらしいですよ!買うしか無い!買うしか無い!(回し者)

アニメには登場しないキャラなので(狐の姿なら毎回出てるけど)、当然衣装とかは売ってません。というわけで、人生初のコスプレは衣装から全部自作です!何故人はマイナーなキャラクターのコスプレをしたくなってしまうのか!!

洋裁をするのは小学校以来で、その時は間違えて縫って履けないズボンを作って、そのまま心折れて提出もせずに卒業したような気がします。そんな感じです。あと3年くらい前に女装してニコ生でiPadケース縫った気がしますが、ありゃノーカン。

作る時に考えたことを、断片的にメモっておこうと思います。

服の構成要素は既成品改造レベル

コンちゃん(人型)の可愛さに気づいて現在パソコンの背景にしている画像がこちらです。

20140520_4この画像をよく見ると、全体の衣装はファンタジーな感じでとてもかわいい1のですが、個々の要素は既存の服をすこし変えたものであることが分かります。

  • マントとフードは、レインポンチョのマントの形を変えたりフードに耳生やしたりすればいけそう
  • 上着は、和服から袖を取ったもの。帯は普通に帯。
  • ズボンは普通の半ズボンの上の部分を袴みたいにすればそれっぽくなりそう。

この時点で、既存の洋裁の型紙を改造すれば、無理をせず十分それらしい衣装が作れそうだと分かります。

あとさらに洋裁初心者には大変うれしいと思ったポイントが

  • とりあえず赤い布と白い布さえあれば作れる

というところです。たくさん材料があると混乱してミスをしてしまいそうなので、とりあえず白い布と赤い布を大量に買っておけばあとは形をどうしてどう縫うかの問題に落とし込めるのは大変好都合に見えました。

 ミシンは適当で良かったっぽい

ミシンすら持ってなかったので、適当に買いました。

ミシンもたくさん種類が出ていてどれ買えばいいの~!?という感じだったのですが、ふなっしーのコスプレをしている超上級者の友人に「3万くらいのコンピュータミシンならどれでもOKじゃない?」みたいなアドバイスをしてもらったので、適当に安いのを選びました。実際使ってみて過不足なかったので、大丈夫だったっぽい。

フットペダルはやっぱり便利ですね。あと最近の機種はみんなそうなのかもしれませんが、この機種は糸通しがミシン本体に付いてて、糸を針に通すのがめっちゃ楽でした!

ちなみにこのミシンが来たのは木曜日で、学園祭は土曜日からだったので非常にタイトなスケジュールでした。もっと計画的になろう(◞‸◟ )

既成品の型紙を集める

大学近く(といっても徒歩20分)にユザワヤが最近できたので、型紙含め材料全部買ってきました。

麻の入った朱色の布がたまたまあったのがすごいラッキーでした。結構余ったので、巫女装束とか作っても良いかな~?

制作

しっぽ⇨袋を作ってフェイクファーをはりつける

しっぽに関しては、四角い袋を適当にぬって先端を丸くしたあと、フェイクファーの帯をボンドで貼り付けまくりました。

SONY DSCそれっぽくなったものの、ボンドが固まっているので、触ると硬くなってしまってモフれないのが難点です。布状になってるフェイクファーがあれば、直接それを袋にすることも可能だったかも…。

マント⇨ほぼそのままだけど裏地大変

マントは今回の衣装にいらなさそうなディテールを適当に削りつつ、マントの形を型紙レベルで書き換えただけです。一番大変だったのは赤い裏地を付ける所で、そのために赤い布と白い布の両方で同じ型紙をカットして縫い合わせないといけないので、作業が普通に作る時の倍以上になります。

SONY DSC装飾とかは全部ボンドで貼り付けました。ボンド便利!!上の写真だと頭にちゃんと葉っぱがついてますが、この葉っぱはリアル葉っぱです(つくりものの葉っぱを買い忘れた)

ズボン⇨袴風にするためにサイドをカットする

SONY DSC袴は昔、高校(ほとんど行ってない)で剣道の授業があって着たことがあったので、それを思い出して、適当に腰のところをカットして、ゴムじゃなくて帯で腰を縛るようにしています。本当の袴はさらに前のほうは5本のひだが付いていたりするのですが、コンちゃんの衣装はそうではなさそう。

結局上着のおはしょりで隠れるのであまり意味はありませんが、しっぽがズボンに固定されていて重いので、それをしっかり持ち上げるのにはちょうど良かったです。

足が出てくるところ(ボタンがついてるところ)がすこしくびれているのがとても可愛いので、再現しようと先の方だけ細くなるように、足先につれて太くなっていた所を途中から細くなるように型紙を改造したのですが、あまりそれっぽい感じにはなってくれませんでした…。

上着⇨おはしょりを作るために長くした

上着は二部ゆかたの上の部分だけ利用することで制作しました。袖も付けていません。

SONY DSCこの二部ゆかたの型紙では、「おはしょり」は浴衣が上下に分かれていることで表現されているのですが、今回は普通に「おはしょり」を作りたかったので、上着の型紙を30cmくらい延長して、長めにしました。

「縫う」より「切る」方が大事

今回やってみて感じたのは、ミシンを使うからついつい「縫う」のが大事なように思ってしまうけど、それよりもちゃんと綺麗に型紙を布に移して、綺麗にカットするのが大事だったということです。華々しさがない作業なのでついつい手を抜きたくなってしまうのですが、ここで手を抜くと後の作業全てに響きます。

小学生の時の私はまち針とかもめんどくさくて殆ど使わなかったような記憶があるのですが、これもちゃんとまち針さして布を固定するのは非常に大事でした。

反省点

反省点もいくつか有ります。

  • 子供用のレインポンチョの型紙を使ったせいで、頭とマントの接合部が首より短くなって着れなくなってしまいました。
    無理やり布でつぎたしましたが、正直不格好。
  • 浴衣の構造を理解せず作っていたので、襟の縫い代が外に出ています。まさか、襟が胸のところまで続いてるなんて知らなかったんだもん!
  • ズボンは長すぎたと思います。足が出てくるところのくびれもとてもかわいいと思うのですが、ここがうまく再現できなかったのは残念です。

あと中の上着に袖が全くありませんが、本当はあります。

20140520_5コミケもこのコスプレしようと思ってますが、その時までに改善できるのか…!?

あと、ちさちゃん(1/6人形)と同じ格好をしてみるのが夢(?)なので、1/6サイズとかで作ってみてもいいかも!

とりあえずTwitterで「作ってみようかな」とつぶやいてみる効果

今まで、「作る」って周りに宣言しちゃうと作った気になるからあんまり言わないようにしてたのですが、今回はこの一言で「ぜったい作るぞ~」みたいな気分になってきたので、「作ろうかな」という思いつきを具体的な言葉にしてつぶやきまくってみるのは結構有用かもな~と思いました。あとは一日24時間の制約と肉体的疲労にどう打ち勝つか…(◞‸◟ )

誰も知らないキャラのコスプレってどういう行為なんだ

コスプレって、例えば初音ミクとか有名なキャラの格好をして「あっ初音ミクだ!似てる!」って思ってもらう一種の相互的なコミュニケーションだと思うんですけど、コンちゃんはほとんど誰も知らないので(とはいえ聞かれたので説明したら2人くらい分かってもらえた)、コミュニケーションが成立してない気がするんですけど、

私が満足してるのでそれでいいとおもいます(小学生並の感想)

  1. 特にマントの裏地が赤いのがほんとーにかわいーとおもいます! []

32bit GCCのunordered_mapに下位32ビットが同じキーを送り続けると死ぬ

Posted on

前回の記事の最後で「unsigned long long intのハッシュは(size_t = 32ビットの環境で)下位32ビットしか見てない」みたいな事を書きました。

これが本当にそうだとすると、この挙動に依存しているunordered_mapに、下位32ビットが同じ値のキーを送り続けるとアルゴリズムの教科書通りの最悪ケースになって悲惨なことになるのでは…!という事が容易に予想されます。

ということで、やってみました。

実験環境

  • Windows 7 64bit
  • MinGW 32 g++ 4.8.1

ソースコードはこちら。

#include <unordered_map>
#include <cstdio>
#include <chrono>

unsigned long getT()
{
    return std::chrono::system_clock::now().time_since_epoch() / std::chrono::milliseconds(1);
}

int main()
{
    {
        std::unordered_map<unsigned long long int, int> map;
        unsigned long t = getT();
        for(int j=0;j<100;++j)
        for(unsigned long long int i=0;i<0x1000;++i){
            map.insert(std::make_pair(i << 17, 10));
        }
        printf("Case 1 Time: %d\n", static_cast<int>(getT() - t));
		fflush(stdout);
    }
    {
        std::unordered_map<unsigned long long int, int> map;
        unsigned long t = getT();
        for(int j=0;j<100;++j)
        for(unsigned long long int i=0;i<0x1000;++i){
            map.insert(std::make_pair(i << 33, 10));
        }
        printf("Case 2 Time: %d\n", static_cast<int>(getT() - t));
		fflush(stdout);
    }
}

シフトする値だけ少し変更することで、他の要素の影響を出来る限り減らしました。

結果

g++ -o test -std=gnu++11 test.cpp
./test
Case 1 Time: 120
Case 2 Time: 67825

というわけで500倍くらい遅かったので、本当に線形探索してるっぽいです。本当にありがとうごいました。

64bit環境なら問題ありません

64bit環境のwandboxで実行すると、64bit環境ではsize_tも64ビットになってくれるので、

Start

Case 1 Time: 287
Case 2 Time: 279

0

Finish

というわけでほぼ同じ時間で終わります。


gccのunordered_mapの実装を読んでみる

Posted on

最近諸事情でアルゴリズムイントロダクションを読んでおります。

で、平均挿入・検索・削除時間がO(1)である謎のテクノロジー、ハッシュテーブルの章を読んだので、せっかくなのでSTLの実装を調べてみました。

なお、コンパイル時の条件によって色々実装が分岐するようですが、面倒臭いので適当に目についたコードを読んでいます。なので、組み合わせとして一緒に呼ばれないコードを読んでいる可能性がありますが、ご了承下さい。

とりあえず、挿入を足がかりに

std::unorderd_mapの実装は、いろいろたらい回しにされたあとbits/hashtable.hに記述されているらしいことが分かります。とりあえず挿入っぽい処理を調べてみると、こんな感じになります。

  template
    void
    _Hashtable<_Key, _Value, _Alloc, _ExtractKey, _Equal, 	       _H1, _H2, _Hash, _RehashPolicy, _Traits>::
    _M_insert_bucket_begin(size_type __bkt, __node_type* __node)
    {
      if (_M_buckets[__bkt])
	{
	  // Bucket is not empty, we just need to insert the new node
	  // after the bucket before begin.
	  __node->_M_nxt = _M_buckets[__bkt]->_M_nxt;
	  _M_buckets[__bkt]->_M_nxt = __node;
	}
      else
	{
	  // The bucket is empty, the new node is inserted at the
	  // beginning of the singly-linked list and the bucket will
	  // contain _M_before_begin pointer.
	  __node->_M_nxt = _M_before_begin._M_nxt;
	  _M_before_begin._M_nxt = __node;
	  if (__node->_M_nxt)
	    // We must update former begin bucket that is pointing to
	    // _M_before_begin.
	    _M_buckets[_M_bucket_index(__node->_M_next())] = __node;
	  _M_buckets[__bkt] = &_M_before_begin;
	}
    }

どうやら、衝突の解決には片方向リストを使った方法を使っているようです。アルゴリズムイントロダクションには「削除が必要な時は片方向リストを使った方法がいい」みたいなことが書いてあったような気がするのですが、そのとおりっぽいです。あと番兵は使ってないみたいですね。やっぱりメモリの節約が優先なのでしょうか。

バケットのインデックスは?

で、この__bktはどこから来てるのかを調べてみると、_M_bucket_indexという関数の値が実際には使われており…bits/hashtable_policy.hに次のようなコードがあります。(型の条件によって複数あるようですがこれが一番なんか読みやすそうなのでこれにしました)

      std::size_t
      _M_bucket_index(const __node_type* __p, std::size_t __n) const
	noexcept( noexcept(declval()(declval()))
		  && noexcept(declval()((__hash_code)0,
						    (std::size_t)0)) )
      { return _M_h2()(_M_h1()(_M_extract()(__p->_M_v())), __n); }

というわけで、h1に値を突っ込んでハッシュを得た後、さらにh2にそのハッシュとバケットの総数を入れてハッシュをバケット数以下にして、実際のバケットのインデックスを求めているみたいです。

__nはバケット数の最大値です。_Key&は例えばstd::unordered_map<string, int>ならstringの方です。_M_h1と_M_h2は実はbits/unorderd_map.hの最初の方にすでに書いてあって、

  template,
	   typename _Pred = std::equal_to<_Key>,
	   typename _Alloc = std::allocator<std::pair >,
	   typename _Tr = __umap_traits<__cache_default<_Key, _Hash>::value>>
    using __umap_hashtable = _Hashtable<_Key, std::pair,
                                        _Alloc, __detail::_Select1st,
				        _Pred, _Hash,
				        __detail::_Mod_range_hashing,
				        __detail::_Default_ranged_hash,
				        __detail::_Prime_rehash_policy, _Tr>;

h1の実体はstd::hash<_Key>(_key const& key)で、これは後で調べますが、とりあえず色々な型についてhash値を返してくれるSTLの公開の標準ライブラリです(皆様もお使いになれます!)

h2は__detail::_Mod_range_hashingで固定で、ユーザーが指定したりは出来ないようです。

h2はただのmod

_Mod_Range_hashingはbits/hashtable_policy.hに実体があり、

  // Many of class template _Hashtable's template parameters are policy
  // classes.  These are defaults for the policies.

  /// Default range hashing function: use division to fold a large number
  /// into the range [0, N).
  struct _Mod_range_hashing
  {
    typedef std::size_t first_argument_type;
    typedef std::size_t second_argument_type;
    typedef std::size_t result_type;

    result_type
    operator()(first_argument_type __num,
	       second_argument_type __den) const noexcept
    { return __num % __den; }
  };

というわけで、散々たらい回しにされた挙句ただのmodでした。万能ハッシュとか使ってるのかと思ったのですがそんなことはなく、h1の実体であるstd::hashに衝突の回避とかの仕事はほぼ丸投げして、バケット数とmodをとっているという感じになります。

std::hashの実体

std::hashの実体は色々なところに散らばっているのですが、std::stringの実装がやっぱり一番気になるのでそれを調べてみると、basic_string.hに実体があって、

  /// std::hash specialization for string.
  template<>
    struct hash<string>
    : public __hash_base<size_t, string>
    {
      size_t
      operator()(const string& __s) const noexcept
      { return std::_Hash_impl::hash(__s.data(), __s.length()); }
    };

  template<>
    struct __is_fast_hash<hash<string>> : std::false_type
    { };

となっております。ほかのstd::wstringとかutf16とかも全部同じ実装です。hashの先を見れば分かるのですが、この関数はvoid*を取る関数なので、文字列としての特性はとくに使わず、バイト列として解釈してhashを計算しているようです。

さて、_Hash_impl::hashはfunctional_hash.hの関数です。このseedの値の意味がよくわかりません…。

    static size_t
    hash(const void* __ptr, size_t __clength,
	 size_t __seed = static_cast<size_t>(0xc70f6907UL))
    { return _Hash_bytes(__ptr, __clength, __seed); }

なんとこの先はlibsupc++という別ライブラリに投げれていて、

  inline std::size_t
  unaligned_load(const char* p)
  {
    std::size_t result;
    __builtin_memcpy(&result, p, sizeof(result));
    return result;
  }

}

namespace std
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION

#if __SIZEOF_SIZE_T__ == 4

  // Implementation of Murmur hash for 32-bit size_t.
  size_t
  _Hash_bytes(const void* ptr, size_t len, size_t seed)
  {
    const size_t m = 0x5bd1e995;
    size_t hash = seed ^ len;
    const char* buf = static_cast<const char*>(ptr);

    // Mix 4 bytes at a time into the hash.
    while(len >= 4)
      {
	size_t k = unaligned_load(buf);
	k *= m;
	k ^= k >> 24;
	k *= m;
	hash *= m;
	hash ^= k;
	buf += 4;
	len -= 4;
      }

    // Handle the last few bytes of the input array.
    switch(len)
      {
      case 3:
	hash ^= static_cast<unsigned char>(buf[2]) << 16;
      case 2:
	hash ^= static_cast<unsigned char>(buf[1]) << 8;
      case 1:
	hash ^= static_cast<unsigned char>(buf[0]);
	hash *= m;
      };

    // Do a few final mixes of the hash.
    hash ^= hash >> 13;
    hash *= m;
    hash ^= hash >> 15;
    return hash;
  }

この関数もなんとなくバイト列中のたくさんのビットを反映させていることはなんとなく分かるのですが、謎のマジックワードもあるし、謎です。

やっぱり本と実装はすこし違いますね~とおもったのでした。おしまい。

他の型のハッシュは?

同じfunctional_hash.hにマクロがあって、

  // Explicit specializations for integer types.
#define _Cxx_hashtable_define_trivial_hash(_Tp) 	\
  template<>						\
    struct hash<_Tp> : public __hash_base<size_t, _Tp>  \
    {                                                   \
      size_t                                            \
      operator()(_Tp __val) const noexcept              \
      { return static_cast<size_t>(__val); }            \
    };

  /// Explicit specialization for bool.
  _Cxx_hashtable_define_trivial_hash(bool)

  /// Explicit specialization for unsigned long long.
  _Cxx_hashtable_define_trivial_hash(unsigned long long)

という感じです。なんとstatic_castするだけ!型に収まってるならともかく、unsigned long longまでstatic_castしちゃうのは本当にそれでいいの?感があります。…ということは32bit環境下でunordered_mapに下32ビットが同じunsigned long longを送り続けると死ぬほど遅い可能性がある…? →実際に実験したら遅かった

ちなみにこのコードのすぐ下にありますが、double/floatに関してはバイト列として文字列と同じハッシュ関数を使っています。うーん…(^_^;)


bashで3Dプログラミング

Posted on

今日はbashと基本的なシェルコマンドだけで3Dレンダリング(ワイヤフレーム)をしてみました!


Github: ledyba / Bash3D

bashプログラミングテクニック

今回つかったbashのプログラミングテクニックについて、記録のついでに少しまとめます。シェルスクリプトの基本的な所は書きません。

あと、3Dプログラミングに関しては本に書いてある事以上のことが書けそうにないので、(今は)書きません。

bashのデータ構造

bashでそこそこ本格的にプログラミングする事を考えた上でまっさきに「辛そう…」ってなるのはデータ構造ですよね!今回のソフトでは、ベクトルや行列を綺麗に扱えないと困ってしまうでしょう。まずはそこを見ます。

bashの「値」は基本的に全てただの文字列ですが、「リスト構造」にだけはネイティブで対応しています。こんな感じ。

$ lst=(1 2 3)
$ echo ${lst[0]}
1
$ echo ${lst[1]}
2
$ echo ${lst[2]}
3
$ echo ${lst[3]}
#空行

データにアクセスする時に、必ず${}でくくらないといけないのが少し面倒です。ちなみに添字を付けないとリスト全体…と思いたいところですが、最初の要素が出てきます。

$ echo $lst
1

全体を表示するには、${lst[@]}とします。

$ echo ${lst[@]}
1 2 3

リストを定義するときに使ったカッコが出てこないので、リストをコピーしたい場合はそのカッコを補わないといけません。

$ lst2=(${lst[@]})
$ echo ${lst2[@]}
1 2 3

慣れればなんとかなるレベルですが、面倒です。

さらに気を付けないといけないのは、”(1 2 3)”という文字列を返す関数を$(command)で受けて変数に代入してもリストになりません。

$ function fun() { echo \(1 2 3\); }
$ lst=$(fun)
$ echo $lst
(1 2 3)

ただの文字列として代入されます。どうしてもしたい場合は、evalを使わないといけません。

#文字列展開が起こってから実行されるので、lst=(1 2 3)というコマンドが実行されることになる
$ eval "lst=$(fun)"
$ echo ${lst[1]}
2
$ echo ${lst[@]}
1 2 3

通常の用途ならそれで済むのですが、今回は出来る限り高速に処理したいので、こうやってコマンドが沢山必要になるのは避けたいところです。

ところで。

後に言うダイナミックスコープとevalを組み合わせると連想配列みたいなことも出来るのですが、今回のプログラムでは基本的にこのリストでデータ構造を管理しています。

インデックスを指定したアクセスをするためにコストが高いevalを行わなくてよいこと、ベクトルと行列が主役なのでリストで十分な事が理由です。

$ matrix=(MAT 1 0 0 0 0 1 .....)
$ vector=(VEC 0 0 0 1)

最初の要素に型の名前を入れて適宜assertすることで、動的型検査みたいなことをしています。

変数のスコープ

名前空間は基本的にひとつしかありません。関数のネストは一切関係なく上書きされていきます。

$ val=0
$ function fun() { val=1; }
$ fun
$ echo $val
1 #関数内の代入だが、グローバル領域が上書きされてしまった
$ function fun2() { val2=2;echo $val2; }
$ fun2
2
$ echo $val2
2 #オフサイドルールではなく、存在しない変数はグローバルに作られるみたい

一種のダイナミックスコープみたいなものだと理解しておけば大体書けます。

呼び出し元の変数がうっかり書き換えられてしまう可能性があるのですごく不便ですが、逆に利用することも出来ます。

「関数の引数に変数名を渡すと、evalと組み合わせることで指定した変数に結果を格納する関数」が作れます。こんな感じ:

$function setOne() { eval "$1=1"; }
$ setOne val1
$ echo $val1
1
$ setOne val2
$ echo $val2
1

これを使って「最初の引数は入力、最後の引数は出力」と規約を決めて関数を設計しています。たとえば行列の計算はこちらにあります。見た目がアセンブラみたいになりますネ。

「こんなスコープ再帰が出来ないのでは?」という気がしますが、案外なんとかなります。というのも、$(command)はサブシェルとかいうので実行されるらしくて、さらにサブシェルは別の変数空間で動作するようです。つまり、このスコープの問題は起こりません。…たぶん。よく分からない(◞‸◟ )

 #再帰関数といえば階乗である
function kaijo() {
    arg=$1
    next=$(expr $arg - 1)
    if [ $arg -gt 1 ]; then
        result=`kaijo $next`
        echo $(expr $arg \* $result)
    else
        echo 1
    fi
}

kaijo 10 --> 3628800

ついでにevalとこのスコープの仕様を使うと、擬似的な連想配列を作れますが、値のコピーが難しいこともあって使いませんでした。

$ val_key1=123
$ val_key2=456
$ function showKey1() { echo $(eval "echo \$$1_key1"); };
$ showKey1 val
123

固定小数点とexprコマンド

シェルスクリプトで計算をするには、exprコマンドを使う必要があります。このコマンドは整数の計算にしか対応していません。

$ expr 1 \* 2 \* 3
6
$ expr 1 \* 2 \* 3.2
expr: non-numeric argument

3Dソフトでそれでは困るので、整数だけで小数を計算します。それにはどうするかというと、扱う小数を適当に1000倍とかして、それを小数に見立てます。つまり、1000は1、1001は1.001、3140は3.140だと思うことにします。その場合、計算する時にちょっとした手間が必要になります。

#3.14 + 3.14
$ expr 3140 + 3140 #足し算・引き算はそのまま足し算で大丈夫です。
6280
#3.14 / 3.14
$ expr 3140 \* 1000 / 3140 #割り算するときは、倍率で掛けてから割ります。この世界では「1.0」は1000であることに注意
1000
#3.14 * 3.14
$ expr 3140 \* 3140 / 1000 #掛け算するときは、掛け算してから倍率で割ります。
9859

倍率を1000じゃなくてもっと高くすれば、もっと小数演算の結果が正確になります。

私の環境だけかもしれませんが、exprはどんな大きな数でもちゃんと扱ってくれるようなので、ガンガン倍率を上げていけば誤差がまったく気にならなくなります。これはすごい楽…(楽すぎた)

なお、sin/cosは(テイラー展開を自分で展開しない限り)計算できないようだったので、0~360度について事前に求めておいて、これを先ほどのリスト構造にいれて使っています。平方根も使いたかったんですが、どうやっても無理そうなので諦めました。

表示:オフスクリーンレンダリング

もちろんターミナルにはピクセルがないので、ターミナル上に表示する文字で代用します。普通に表示するなら、tput cup (line) (col)を使うとターミナル上の任意の位置に移動できるので、これを使いたいところです。

# 斜め45度に線を引いてみた
for ((x=0;x<10;++x));do
    for ((y=0;y<10;++y));do
        tput cup y x
        echo -n \*
    done
done

が、この方法では線を引く様子が目で分かるぐらい時間がかかります。これで最初レンダリングしたみたのですが、レンダリング中に線を引いている様子がバッチリ見えてしまい、あんまり回転アニメーションしてるように見えませんでした。

ほんもののゲームでもそういう、更新途中の絵が表示されてしまってちらつく問題があります。ほんもののゲームでは、そういう時は、書いてもすぐには表示されない「裏の画面」を用意しておいて、そちらに時間を掛けてレンダリングして、最後に表示される「表の画面」と入れ替える「オフスクリーンレンダリング」と「ダブルバッファリング」をしているのでした。私も同様にそのテクニックを使いました。

具体的には、こんな感じ。

SCR_BUFF=$(head -c $size < /dev/zero | tr '\0' '_')

今回の「ピクセル」は文字ですから、「裏の画面」はただの文字列になります。
画面のサイズを取得しておいて、そのサイズ分何もない文字列を意味する「_」をtrを使って入れておきます。空白でいいじゃないかと思うかもしれませんが、シェルスクリプトでは空白はかなり重要な意味を持っていて面倒なのでこうしています。

「ドットを打つ」時は、対応する文字列中のインデックスの文字を*に書き換えます。

SCR_BUFF=$(echo $SCR_BUFF | sed s/./*/<置き換える文字の位置>)

zshだとstr[12]=\*ってやると文字列を更新できるのですが、bashでは使えなかったので、sedを使って書き換えています。この処理が非常に重く、今回のレンダリングエンジン(?)の処理時間の殆どはこのドットを打つ処理に費やされています。上手い方法を知ってたら教えてください!

最後に表の画面に更新する方法ですが、単純にechoするだけです。この時、何も入っていない所には_が入っていたので、空白に直します。

$ tput cup 0 0
$ echo -n $SCR_BUFF | tr '_' ' '

ちゃんと幅に合わせて改行されるので、うまく表示できます。

高速化

インライン展開 → プロセスforkを抑える

当初綺麗に関数を使って書いた所、行列の掛け算に1秒近く掛かってしまい、さすがにレンダリングは出来ないのかな…と諦めかけたのですが、行列やベクトルの掛け算や足し算を全てインライン展開したところ、座標計算に関しては十二分すぎる程の速度が出せるようになりました。前述した通り、今回遅いのはオフスクリーンにピクセル(文字)を書き込む処理が異常に時間が掛かっているからです。

やはりプロセスの起動はものすごい重い処理のようで、可読性をそこまで損なわない範囲でプロセスの呼び出しを一箇所にまとめています。

感想

線を引くのにえげつない時間が掛かるのを見た時、初めて触ったBASICの処理系を思い出しました。懐かしい。

ところで、なんで3Dライブラリを作るとOpenGL風の名前にしちゃうんだろ。

ニコ動に「次はシェーディングだな!」というコメントがついているのですが、テクスチャの表示をやってみたいです!7色くらいしか使えないので、たいへんそう…。影は「暗い色」が基本的に無いのでかなりつらい気がします。


東北地方太平洋沖地震の地震波を音にしてみた

Posted on

地球物理をやる学科に居るので、地震も学んでいるのですが、ふと思いついたので地震計の観測した波形を音に見立ててWAVファイルにしてみました。

気象庁の強震波形(平成23年(2011年)東北地方太平洋沖地震)のページから、福島県郡山市朝日の、震度6のデータです。なんだか、雷みたいな音がしますね!

ソースコードはこちらです。PythonでCSV読んでWAVに書いているだけですが、一応使い方のサンプルとしては使えるのではないでしょうか。


凪のあすからの恋愛関係からページランクを求めてみました

Posted on

最近友人に薦められて、今日が最終回らしい「凪のあすから」を見ています。今6話なのですが、恋愛関係が大変複雑です。「深夜の昼ドラ? アニメ『凪のあすから』の人物相関図つくったよ」というページがあったので、引用しますと、こんな感じ。

うーん、グラフ理論って感じだ。そう思ってからピクシブの非公式カップリングのページを見てしまうと、もはや隣接行列か何かにしか見えなくなってきた…。

20130403隣接行列…グラフ理論…はっ!ページランク求めようページランク!

というわけで、今まで一回ぐらいやってみようと思っていた、ページランクの計算を実際にやってみました。こちらの「Google の秘密 – PageRank 徹底解説」を参考にしました。人間をページに見立てて、好意をリンクに見立てます。

ページランクとは結局何なのか

昔私が読んだ本だと、

  • 全てのページはページランク値PRを持っていて、
  • 他のリンクしているページそれぞれににPR/N(Nはページ内のリンク総数)を与えており、
  • ページランクの値PRはリンクされた他のページから貰ったページランク値の和に等しい

というのを読んだ気がします。つまり、貰ったページランク=あげたページランク=そのページのページランクという等式が成り立つということです。

元の論文にもそんな図が載っています。

20130403-2

Lawrence Page, Sergey Brin, Rajeev Motwani, Terry Winograd, ‘The PageRank Citation Ranking: Bringing Order to the Web’, 1998,
http://www-db.stanford.edu/~backrub/pageranksub.ps

これをグラフ理論と確率過程の話で表現すると、こうなります。

  • 各ページから他のページへ移動する確率の表を用意する(とりあえずすべてのリンクで等確率=1/リンク総数とする)
  • その表に従ってあるページから出発してリンクをたどり続ける
  • すると、どのページから出発しても、N回リンクをたどった時にあるページPに居る確率がN+1回目のそれと変わらなくなってくる(これを定常状態と呼びます)
  • その時の、N回目でそれぞれのページに居る確率がページランクです。

これがホントに元の話と一緒なのか?ということなのですが、

  • “もらうページランク”=(自分も含めた他のページにN−1回目に居る確率*それぞれの遷移確率)の和=N回目にそのページにいる確率
  • “あげるページランク”=(N回目にそのページにいる確率 * そのページから他のページへの遷移確率)の和=N回目にそのページにいる確率1

なので、もらうページランク=あげるページランクという関係は確かに成り立ちます。

求めてみる

ページ遷移確率表=行列

実際にもとめてみましょう。まずは、鍵となるページ遷移確率表を用意します。さっきの表を見ながら、ひーくんからまなかへ遷移する確率は1…みたいな感じで行列を作るとこうなります。一応、今回は友情は無視して恋愛関係だけ考えます。

M =
 ひーくん ちさき  まなか   要    さゆ   美海   つむぐ  あかり  至    みをり ←好き ↓好かれる
   0.00000   1.00000   0.00000   0.00000   0.00000   1.00000   0.10000   0.00000   0.00000   0.00000 ひーくん
   0.00000   0.00000   0.00000   1.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 ちさき
   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 まなか
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000   0.10000   0.00000   0.00000   0.00000 要
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 さゆ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 美海
   0.00000   0.00000   1.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 つむぐ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 あかり
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   1.00000   0.00000   1.00000 至
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 みをり

普通、ページには複数のリンクがありますが、このアニメだと思い人が二人以上居るのは不倫と勘違いされて殴られた至だけです。

あと6話時点ではつむぐは誰も好きではないので、つむぐの行が全部0になるはずなのですが、そうしてしまうとページランクが求まりません2

これを元論文ではdangling pageと呼んでいます。元論文では「あんまり重要じゃないから削除して計算しています」と無視してしまっているのですが、つむぐくんを無視するのはかわいそう(?)なのと、仮につむぐくんを消してしまうと今度はまなかがdangling pageになって…となって最終的に大人組しか居なくなってしまうので、次に紹介するrandom surferモデルに習って、みんなを平等に好きなことにしました。

計算は掛け算するだけ

元の定義通り、「あるページからスタートしてリンクをたどり続けた時に、それぞれのページに居る確率が一定になってくる」まで遷移させ続けます。

たとえば、ひーくんからはじめましょう。はじめはひーくんに居ると仮定しているので、ひーくんに居る確率は1で、他みんな確率0とします。

1回リンクをたどると、グラフを見るとひーくんはまなかが好きなので、まなかの確率が1になるはず。

これを我々が先ほど作った行列を使って計算しようとするなら、先ほどの行列Mに、ひーくんを1としたベクトルPを掛ければよいです。

octaveでやるとこんな感じです。

octave:95> M=[ 0 1 0 0 0 1 1/10 0 0 0 ; 0 0 0 1 0 0 1/10 0 0 0; 1 0 0 0 0 0 1/10 0 0 0; 0 0 0 0 1 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0; 0 0 1 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 1/2 0; 0 0 0 0 0 0 1/10 1 0 1; 0 0 0 0 0 0 1/10 0 1/2 0 ]
(略)
octave:96> P = [1; 0; 0; 0; 0; 0; 0; 0; 0; 0] ←ひーくんに居る確率を1としたベクトル
octave:97> M*P
ans =

   0
   0
   1
   0
   0
   0
   0
   0
   0
   0

というわけで、まなかは三番目の要素で、それが1となったので、見事正しい結果を得られています。これを沢山続けてみましょう。10000回ぐらい掛けると変わらなくなってきます。

octave:98> (M^10000)*P
ans =

   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.22222
   0.55556
   0.22222

…あれぇ?最後の3人、つまり大人組以外ゼロ?なんか間違ったかなぁ?

ランダム・サーファー・モデル

こうなってしまった原因は、大人組から子供組へ繋がっていないことです。つむぐくんは全員が好きだと仮定したので、子供組から大人組へはつむぐくんから遷移することが出来るのですが、その後は大人組の間でぐるぐるぐるぐる巡り続けることになります。ひーくんからはじめても長い目で見ればいつかは大人組に入ってしまい、そこで定常状態となるので、子供組の確率がすべて0になってしまったのです。グラフ理論の言葉でいうと、グラフを辿ると元の位置に必ず戻ってこれる大人組のような部分を再帰類と呼び、そうではなく元に戻れない子供組のような部分を非再帰と呼ぶそうですが、あまりこれ自信をもって解説できないので、しません…。

元の論文ではこれをRank sinkと読んでいて、そのためにランダム・サーファー・モデルというのを考慮しています。これはどういうモデルかというと、「時々はまったくリンクされていない別のページに飛ぶ」と考慮しているモデルです。具体的には、新しい遷移確率行列を次のように定義します。

octave:99> M2=(M*1/1.15)+((1/10)*0.15/1.15)

これは、まず今までの確率を少し下げた上で、あるページから他のページへ僅かな確率で飛ぶようにした、ということです。凪のあすからの言葉を使えば、ひょっとするとひーくんは他の女の子(や男の子)にも興味を持つかもしれないよねというモデルです。

この修正したモデルでは、あるページから他のページへ少ない確率ながら必ずジャンプできるので、Rank sinkのような問題は起こりません。

早速、計算してみましょう。

octave:100> (M2^10000)*P
   0.095976 ひーくん
   0.060684 ちさき
   0.106569 まなか
   0.043208 かなめ
   0.023111 さゆ
   0.023111 美海
   0.115780 つむぐ
   0.135981 あかり
   0.259599 至
   0.135981 みをり

やはり両思いは強いのですね!(?) 子供組を見ると、やはり特に重要そうなキャラクターであるひーくん、まなか、つむぐの三人のページランクが高いので、おおよそページランクの妥当性を示せたと言っていいんでしょうか!?とりあえず、この三人に注目して今後は視聴を続けたいと思います!

蛇足:固有値・固有ベクトルを使って求める

なんか10000回も行列を掛け算するのは…という感じしますか?元の論文だとその方法で計算してると書いてあるのですが、これは行列の固有値・固有ベクトルを使っても求めることが出来ます。

固有値・固有ベクトルをおさらいしておくと、

  • ある行列Mに対する、
  • 固有値μ・固有ベクトルVは、
  • 次のような関係を満たすものです
  • MV = μV

ここで、固有値μ=1なら、MV=Vとなって、N回目とN+1回目でのそれぞれのページに居る確率が全く同じという事になり、ちょうど私達が求めたかったベクトルと一致しますね!

octaveで求めるとこんな感じになります。

octave:142> eig(M2)
ans =
   1.00000 + 0.00000i ←求めたい、固有値1の固有ベクトルは1番目だと分かる
   0.78950 + 0.00000i
   0.19713 + 0.55213i
   0.19713 - 0.55213i
  -0.86957 + 0.00000i
  -0.59237 + 0.00000i
  -0.25222 + 0.45314i
  -0.25222 - 0.45314i
   0.00000 + 0.00000i
   0.00000 + 0.00000i

octave:143> [V, D] = eig(M2); ←固有値と固有ベクトルを再度格納
octave:144> V(:,1) ←固有ベクトル1番目を取得する
ans =
   0.252078
   0.159383
   0.279899
   0.113484
   0.060701
   0.060701
   0.304091
   0.357146
   0.681825
   0.357146
↑固有ベクトルは整数倍しても固有ベクトルで、octaveはベクトルの大きさを1に正規化するので、さっきの結果と合わない
octave:146> V(:,1) ./ sum(V(:,1)) ←ベクトルの要素の合計が1となるようにすると…
ans =
   0.095976
   0.060684
   0.106569
   0.043208
   0.023111
   0.023111
   0.115780
   0.135981
   0.259599
   0.135981
↑先ほどの結果と一致する!

こんな感じで計算できます。固有値ってこんなのにも使えるんですね〜。ちなみに、この行列の固有値のうち、絶対値が一番大きいものはかならず1になることがペロン=フロベニウスの定理から言えるらしいのですが、よくわからないので解説できません…。

感想

凪のあすからは線形代数と確率過程とグラフ理論を学べる神アニメ

ていうか実は前期のアニメで一番好きだったのは「いなり、こんこん、恋いろは」で、京都と出雲にまで行ってきたのですが、それを書く日は来るのだろうか…!?

  1. 他のページへの遷移確率の和は1です。リンク総数*(1/リンク総数)=1。 []
  2. つむぐから他のページに遷移できないので、確率がそこで皆0になってしまう []