デモはこのページの一番下にある。
ソースコードは yuki67/MinCamlOnline にある。
何をしたか
前の記事 で bucklescript を使って OCaml プログラムをブラウザ上で動かすというのをやった。
この記事で試したのはかなり小さいプログラムだったので、もう少し大きいプログラムでも上手く行くのかどうか気になり試した。
もう少し大きいプログラムとして MinCaml を使った。 MicCaml なら自分もそれなりに中身を理解していて、プレーンの OCaml で書かれていることを知っているのでちょうどいいと思った。
感想
MinCaml をブラウザ上で動かすために bucklescript を設定するのは思っていたよりも簡単だった。
ただ動かすだけなら ocamlyacc と ocamllex を設定して
C で書かれていた一部のコードを OCaml で書きなおすだけで動くようになった。
動くようにした後にテキストエリアやボタンを設定するのが多少手間だった。
ボタンが押されたときにコンパイルする関数を呼べるようにモジュールの構成を変える必要があったり、
MinCaml にコンパイル結果を文字列として返す関数が無かったのでそれも書いたりした。
これらの手間が生じたのは MinCaml がそのような用途のために書かれていないからであり、OCaml や bucklescript の問題ではない。
変換された javascript プログラムは Parsing, Lexing, Format モジュールを含めて正しく動いた。
ただし、 Format.sprintf が常識的な長さの入力に対して Maximum call stack size exceeded を出して動かないことがあった (デモで mandelbrot をコンパイルするとエラーが出る)。
bucklescript の実装が原因だと思うが、対処する方法もわからないのでそのままにした。
デモ
Program
select program After Type Checking
After K-normalization
After Alpha Conversion
After Optimization
After Closure Conversion
Virtual Assembly Code
After Immediate Value Optimization
After Register Allocation
Final Output
SICP 3.5節 にある無限ストリームを C++ で実装した。
それを使ってフィボナッチ数列や素数、Aitken 加速などを実装できた。
ラムダ式とか auto とかをいい感じに使えた、と思う。
あとこの記事を描くのに noweb を使った文芸的プログラミングを初めて使った。
実装
SICP では無限ストリームは car の値と cdr 部分のストリームを計算する関数の組として実装されている。
たとえばnから始まる整数を並べた無限ストリームは
;; n から始まる整数ストリーム
( define ( integers-starting-from n )
( cons-stream n ( integers-starting-from ( + n 1 )))) となる。
これと同じように C++ でストリームを実装しようとするとすごく面倒臭くなる。
C++ にはラムダ式と auto を使って関数を定義するときには再帰が使えないという制限があるので、上の例と似た関数を
// コンパイルできない
auto integers_starting_from = []( auto n ) { return cons ( n , integers_starting_from ( n + 1 )); };
と定義することができず、型(いまの例なら std::functional<Stream<int>(int)> )を毎回書かなければならない。
流石にこれは面倒臭すぎるので、ストリームを配列っぽく実装した。
T 型のストリームは第n項を計算する関数によって初期化され、
一度計算した値は std::unordered_map<size_t, T> に保存される。
さらに、第n項を計算するときに他の項の値を参照できるように、ストリームを初期化するときに与える
関数には自分自身 (*this) が引数として与えられる。
Code Snippet 1:
Stream の定義
#include <functional>
#include <iostream>
#include <unordered_map>
template < typename T >
struct Stream {
private :
// genrator(self, n) はストリームの第n項を返す
// 第一引数には自分自身 (*this) が与えられるので、第n項を計算するときに他の項を参照することができる
std :: function < T ( Stream & , size_t ) > generator ;
// 計算された値を保存するハッシュ
std :: unordered_map < size_t , T > memo ;
public :
Stream () {}
explicit Stream ( std :: function < T ( Stream & , size_t ) > generator_ )
: generator ( generator_ ) {}
// ストリームの第n項を返す
T & operator []( size_t n ) {
if ( memo . count ( n ) == 1 ) {
// 計算済みなら、その値を返す
return memo [ n ];
} else {
// 計算してないなら、計算を行い、結果を保存してから、結果を返す
memo [ n ] = generator ( * this , n );
return memo [ n ];
}
}
};
以下のプログラムではフィボナッチ数列を計算している。
<< Streamの定義 >>
// フィボナッチ数列
static Stream < size_t > fib ([]( auto & self , auto n ) -> size_t {
if ( n < 2 ) {
return 1 ;
} else {
// 他の項を参照する。メモ化は適切に行われる。
return self [ n - 1 ] + self [ n - 2 ];
}
});
int main () {
for ( size_t i = 0 ; i < 10 ; i ++ ) {
std :: cout << "fib[" << i << "] = " << fib [ i ] << " \n " ;
}
std :: cout << "fib[" << 100 << "] = " << fib [ 100 ] << " \n " ;
return 0 ;
}
fib[0] = 1
fib[1] = 1
fib[2] = 2
fib[3] = 3
fib[4] = 5
fib[5] = 8
fib[6] = 13
fib[7] = 21
fib[8] = 34
fib[9] = 55
fib[100] = 1298777728820984005その他の例:
#include <cassert>
<< Streamの定義 >>
// 整数
static Stream < size_t > integers ([]( auto & , auto n ) { return n ; });
// 素数
static Stream < size_t > primes ([]( auto & self , auto n ) {
static int i = 2 ;
for (;; i ++ ) {
auto flag = true ;
// n番目の素数は、n-1番目までの素数全てで割り切れない最小の自然数
for ( size_t j = 0 ; j < n ; j ++ ) {
flag &= ( i % self [ j ]) != 0 ;
}
if ( flag ) {
return i ;
}
}
});
// 双子素数
static Stream < std :: pair < size_t , size_t >> twin_primes (
[]( auto & , auto n ) -> std :: pair < size_t , size_t > {
// primes[i] と primes[i+1] が双子素数であるようなiを返す
static Stream < size_t > twin_primes_idx ([]( auto & self , auto n ) -> size_t {
if ( n == 0 ) {
return 1 ;
}
auto i = self [ n - 1 ] + 1 ;
for (;; i ++ ) {
if ( primes [ i ] + 2 == primes [ i + 1 ]) {
return i ;
}
}
});
auto younger = primes [ twin_primes_idx [ n ]];
return { younger , younger + 2 };
});
// コラッツ予想の操作をしたときに何回で1になるかを返す
static Stream < size_t > Collatz ([]( auto & self , auto n ) -> size_t {
if ( n == 1 ) {
return 0 ;
} else if (( n % 2 ) == 0 ) {
return self [ n / 2 ] + 1 ;
} else {
return self [ 3 * n + 1 ] + 1 ;
}
});
int main () {
assert ( integers [ 0 ] == 0 );
assert ( integers [ 1 ] == 1 );
assert ( integers [ 2 ] == 2 );
assert ( integers [ 3 ] == 3 );
assert ( integers [ 100 ] == 100 );
assert ( primes [ 0 ] == 2 );
assert ( primes [ 1 ] == 3 );
assert ( primes [ 2 ] == 5 );
assert ( primes [ 3 ] == 7 );
assert ( primes [ 4 ] == 11 );
assert ( primes [ 99 ] == 541 );
assert ( twin_primes [ 0 ]. first == 3 );
assert ( twin_primes [ 1 ]. first == 5 );
assert ( twin_primes [ 2 ]. first == 11 );
assert ( twin_primes [ 3 ]. first == 17 );
assert ( twin_primes [ 99 ]. first == 3821 );
assert ( Collatz [ 1 ] == 0 );
assert ( Collatz [ 11 ] == 14 );
assert ( Collatz [ 27 ] == 111 );
assert ( Collatz [ 871 ] == 178 );
std :: cout << "OK \n " ;
return 0 ;
}
数列の加速
SICP にある Aitken 加速の例も実装できたので解説する。
まず \(\frac{\pi}{4}\) に収束する数列
\begin{align*}
a_{n} = \sum_{i=0}^{n} \frac{(-1)^{n}}{2n+1}
\end{align*}
を計算するストリームを作った。
この数列は
\begin{align*}
a_{n} =
\begin{cases}
1 & (n = 0) \\
と書き直せるので、この数列の第 n 項を計算するストリームを以下のように定義した。
Code Snippet 2:
数列Aの定義
#include <cmath>
// pi / 4 に収束する
static Stream < double > A ([]( auto & self , auto n ) -> double {
if ( n == 0 ) {
return 1 ;
} else {
return self [ n - 1 ] + std :: pow ( - 1 , n ) / ( 2 * n + 1 );
}
});
#include <iostream>
<< Streamの定義 >>
<< 数列 Aの定義 >>
int main () {
for ( size_t i = 0 ; i < 10 ; i ++ ) {
std :: cout << "4*A[" << i << "] = " << 4 * A [ i ] << " \n " ;
}
std :: cout << "4*A[" << 100 << "] = " << 4 * A [ 100 ] << " \n " ;
std :: cout << "4*A[" << 1000 << "] = " << 4 * A [ 1000 ] << " \n " ;
std :: cout << "4*A[" << 10000 << "] = " << 4 * A [ 10000 ] << " \n " ;
return 0 ;
}
4*A[0] = 4
4*A[1] = 2.66667
4*A[2] = 3.46667
4*A[3] = 2.89524
4*A[4] = 3.33968
4*A[5] = 2.97605
4*A[6] = 3.28374
4*A[7] = 3.01707
4*A[8] = 3.25237
4*A[9] = 3.04184
4*A[100] = 3.15149
4*A[1000] = 3.14259
4*A[10000] = 3.14169かなりゆっくりではあるが A が \(\frac{\pi}{4}\) に収束しているのが分かる。
この数列 A に対して Aitken 加速を適用することができる。
数列 \(x_{n}\) に対して Aitken 加速を適用した数列 \(x_{n}^{\prime}\) は
\begin{align*}
x^{\prime}_{n} &= x_{n} - \frac{(\Delta x_{n})^{2}}{\Delta^{2}x_{n}} \\
と定義され、\(x_{n}^{\prime}\) は \(x_{n}\) に比べて収束がとても速くなる(ことがある)。
Aitken 加速をストリームを受け取ってストリームを返す関数として実装した。
Code Snippet 3:
Aitken加速の定義
auto Aitken = []( auto & x ) {
return Stream < double > ([ & x ]( auto & , auto n ) {
auto ddx = ( x [ n + 2 ] - x [ n + 1 ]) - ( x [ n + 1 ] - x [ n ]);
auto dx = x [ n + 1 ] - x [ n ];
if ( ddx == 0 ) {
// 分母が0になるときの例外処理
return x [ n ];
}
return x [ n ] - dx * dx / ddx ;
});
};
Aitken 加速によって収束が速くなることを確認した:
#include <iostream>
<< Streamの定義 >>
<< 数列 Aの定義 >>
<< Aitken加速の定義 >>
int main () {
auto AA = Aitken ( A );
for ( size_t i = 0 ; i < 10 ; i ++ ) {
std :: cout << "4*AA[" << i << "] = " << 4 * AA [ i ] << " \n " ;
}
return 0 ;
}
4*AA[0] = 3.16667
4*AA[1] = 3.13333
4*AA[2] = 3.14524
4*AA[3] = 3.13968
4*AA[4] = 3.14271
4*AA[5] = 3.14088
4*AA[6] = 3.14207
4*AA[7] = 3.14125
4*AA[8] = 3.14184
4*AA[9] = 3.14141第10項で小数以下5桁まで正しく計算できているので、収束が加速されたことが分かる。
Aitken 加速を適用した数列に更に Aitken 加速を適用することができ、さらにその数列を加速することもできる。
この再帰的な Aitken 加速をストリームのストリームを使って実装した。
Code Snippet 4:
再帰的Aitken加速の定義
// 再帰的に Aitken 加速されたストリームのストリームを返す
auto AitkenTableau = []( auto & x ) {
return Stream < Stream < double >> ([ & x ]( auto & self , auto n ) {
if ( n == 0 ) {
return Aitken ( x );
} else {
return Aitken ( self [ n - 1 ]);
}
});
};
auto AitkenRecursive = []( auto & x ) {
auto tableau = AitkenTableau ( x );
// 計算した tableau をクロージャの中に入れる
// tableau は計算の中で更新されるので、mutable が必須
return Stream < double > (
[ tableau ]( auto & , auto n ) mutable { return tableau [ n ][ 0 ]; });
};
#include <iostream>
#include <iomanip>
<< Streamの定義 >>
<< 数列 Aの定義 >>
<< Aitken加速の定義 >>
<< 再帰的 Aitken加速の定義 >>
int main () {
auto AAA = AitkenRecursive ( A );
std :: cout << std :: setprecision ( 15 );
for ( size_t i = 0 ; i < 10 ; i ++ ) {
std :: cout << "4*AAA[" << i << "] = " << 4 * AAA [ i ] << " \n " ;
}
return 0 ;
}
4*AAA[0] = 3.16666666666667
4*AAA[1] = 3.14210526315789
4*AAA[2] = 3.141599357319
4*AAA[3] = 3.14159271403378
4*AAA[4] = 3.14159265397529
4*AAA[5] = 3.14159265359118
4*AAA[6] = 3.14159265358978
4*AAA[7] = 3.1415926535898
4*AAA[8] = 3.14159265358979
4*AAA[9] = 3.14159265358979第9項の値は小数以下14桁まで正しい。
参考
Structure and Interpretation of Computer Programs
jupyter notebook で sympy を使うとセルの実行結果を Latex で表示してくれるので便利だが、
ob-ipython で sympy を使うとセルの実行結果の表示が jupyter notebook ほど便利でない。
#+BEGIN_SRC ipython :session :exports both :results output raw drawer
from sympy import *
init_printing()
x = symbols("x")
eqn = Eq(x**2+3*x-4)
print(latex(eqn))
#+END_SRC
#+RESULTS:
:RESULTS:
x^{2} + 3 x - 4 = 0
:END:from sympy import *
init_printing ()
x = symbols ( "x" )
eqn = Eq ( x ** 2 + 3 * x - 4 )
print ( latex ( eqn )) x2 + 3 x - 4 = 0
幸い org-babel には「src block の実行結果を別の src block に通して、その結果を元の src block の結果にする」
という機能があるので、その機能を使えばこの問題を回避した。
引数付きの src block は emacs lisp だと書きやすいので、それを使って書いた。
sympy の出力にエスケープされたカッコをつけて MathJax のソースコードとして認識されるようにした。
#+name: wrapper-latex
#+begin_src emacs-lisp :var text="" :exports none
(concat "\\(" (replace-regexp-in-string "\n$" "" text) "\\)\n")
#+end_src
#+HEADER: :session :exports both :results output drawer
#+HEADER: :eval no-export :post wrapper-latex(text=*this*)
#+BEGIN_SRC ipython
x = symbols("x")
y = symbols("y", cls=Function)
eqn = Eq(y(x).diff(x), sin(x)*tan(y(x)))
print(latex(eqn))
#+END_SRC
#+RESULTS:
\(\frac{d}{d x} y{\left (x \right )} = \sin{\left (x \right )} \tan{\left (y{\left (x \right )} \right )}\)x = symbols ( "x" )
y = symbols ( "y" , cls = Function )
eqn = Eq ( y ( x ) . diff ( x ), sin ( x ) * tan ( y ( x )))
print ( latex ( eqn )) \(\frac{d}{d x} y{\left (x \right )} = \sin{\left (x \right )} \tan{\left (y{\left (x \right )} \right )}\)
だいぶ良くなった。
次に数式の最初に \displaystyle をいれ、さらに数式全体を special block で囲んで見やすくした。
#+name: wrapper-latex2
#+begin_src emacs-lisp :var text="" :exports none
(concat "#+BEGIN_exa\n\\(\\displaystyle"
(replace-regexp-in-string
"# Out\\[.*\\]:\n: '\\(.*\\)'"
"\\1"
text)
"\\)\n#+END_exa")
#+end_src
#+HEADER: :session :exports both :results value raw
#+HEADER: :eval no-export :post wrapper-latex2(text=*this*)
#+BEGIN_SRC ipython
x = symbols("x")
y = symbols("y", cls=Function)
eqn = Eq(y(x).diff(x), (x-2*y(x)+3)/(2*x+y(x)-4))
latex(eqn)
#+END_SRC
#+RESULTS:
#+BEGIN_exa
\(\displaystyle\\frac{d}{d x} y{\\left (x \\right )} = \\frac{x - 2 y{\\left (x \\right )} + 3}{2 x + y{\\left (x \\right )} - 4}\)
#+END_exax = symbols ( "x" )
y = symbols ( "y" , cls = Function )
eqn = Eq ( y ( x ) . diff ( x ), ( x - 2 * y ( x ) + 3 ) / ( 2 * x + y ( x ) - 4 ))
print ( latex ( eqn ))
\(\displaystyle\frac{d}{d x} y{\left (x \right )} = \frac{x - 2 y{\left (x \right )} + 3}{2 x + y{\left (x \right )} - 4}\)
さらに src block を実行するたびに latex 要素を画像として表示するようにした。
使い勝手が jupyter notebook にかなり近くなった。
( add-hook 'org-babel-after-execute-hook 'org-toggle-latex-fragment ) 参考
Header arguments - Org Babel reference card
PAPI を使うと CPU のハードウェアカウンタを使って様々な値を測定できる。
例えば
実行した命令の数
実行した浮動小数点演算の数
パイプラインがストールした回数
キャッシュヒット/ミス の回数
分岐予測の成功/失敗数
を測定できる。測定した値は FLOPS の測定なりパフォーマンスのチューニングなりに使える。
適切に設定すれば CPU 以外のハードウェアのカウンタも読めるらしいが、使ってないのでよくわからない。
インストール
昔はインストールするのにカーネルにパッチを当てる必要があったりして一苦労だったらしい (検索すると出てくる) が、
最新版の Linux カーネルであればインストールは ./configure; make; make install で済んだ。
利用可能なイベントの確認
papi_avail で利用可能なイベントの一覧を確認できる。コンシューマ向けの CPU だと性能測定用のハードウェアカウンタがしょぼいようなので、インストールしたら必ず確認したほうがいい。
自分の環境では以下のようになった。
$ papi_avail
Available PAPI preset and user defined events plus hardware information.
--------------------------------------------------------------------------------
PAPI version : 5.6.1.0
Operating system : Linux 4.4.0-128-generic
Vendor string and code : GenuineIntel (1, 0x1)
Model string and code : Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz (158, 0x9e)
...
--------------------------------------------------------------------------------
================================================================================
PAPI Preset Events
================================================================================
Name Code Avail Deriv Description (Note)
PAPI_L1_DCM 0x80000000 Yes No Level 1 data cache misses
PAPI_L1_ICM 0x80000001 Yes No Level 1 instruction cache misses
PAPI_L2_DCM 0x80000002 Yes Yes Level 2 data cache misses
PAPI_L2_ICM 0x80000003 Yes No Level 2 instruction cache misses
PAPI_L3_DCM 0x80000004 No No Level 3 data cache misses
PAPI_L3_ICM 0x80000005 No No Level 3 instruction cache misses
...
PAPI_FML_INS 0x80000061 No No Floating point multiply instructions
PAPI_FAD_INS 0x80000062 No No Floating point add instructions
PAPI_FDV_INS 0x80000063 No No Floating point divide instructions
PAPI_FSQ_INS 0x80000064 No No Floating point square root instructions
PAPI_FNV_INS 0x80000065 No No Floating point inverse instructions
PAPI_FP_OPS 0x80000066 No No Floating point operations
PAPI_SP_OPS 0x80000067 Yes Yes Floating point operations; optimized to count scaled single precision vector operations
PAPI_DP_OPS 0x80000068 Yes Yes Floating point operations; optimized to count scaled double precision vector operations
PAPI_VEC_SP 0x80000069 Yes Yes Single precision vector/SIMD instructions
PAPI_VEC_DP 0x8000006a Yes Yes Double precision vector/SIMD instructions
PAPI_REF_CYC 0x8000006b Yes No Reference clock cycles
--------------------------------------------------------------------------------
Of 108 possible events, 59 are available, of which 18 are derived.出力の最終行から、全体の約半分 (59 / 108) のイベントしか測定できないことが分かる。
なお PAPI_FP_OPS という浮動小数点演算の数を数えるための一番基本的なイベントが使えないことから、
PAPI のサンプルコードやテストが上手く実行できなかった。
浮動小数点演算の回数自体は PAPI_DP_OPS というイベントを使って測定できたのでよかったが、
一番最初に実行するサンプルコードが segmantation fault で落ちるので原因の調査にかなり時間を食った。
使う
以下のプログラムは浮動小数点演算の回数、CPU時間を計測し、そこから計算される GFLOPS と合わせて表示させた。
それから実行された命令の数も計測して表示させている。
#include <iostream>
#include <papi.h>
#include <vector>
void print_papi_error ( int error_num ) {
std :: cout << "PAPI error " << error_num << ": " << PAPI_strerror ( error_num )
<< " \n " ;
}
int main () {
int retval ;
auto a = 0.0 ;
// 測定するイベント
std :: vector < int > events = { PAPI_TOT_INS , PAPI_DP_OPS };
// 測定結果を収める場所
std :: vector < long_long > values ( events . size ());
// カウンタをスタートさせる
// events が多すぎる場合はここでエラーが出る
retval = PAPI_start_counters ( events . data (), events . size ());
if ( retval != PAPI_OK ) {
print_papi_error ( retval );
}
auto cpu_time_begin_usec = PAPI_get_virt_usec ();
// 浮動小数点演算を 10^8 回行う
// 前後の命令が依存関係を持つので性能は出ない
for ( auto i = 0 ; i < 1.e8 ; i ++ ) {
a += 1 ;
}
auto cpu_time_end_usec = PAPI_get_virt_usec ();
// カウンタの値を読む
// 値を読んでもカウンタは止まらない
retval = PAPI_read_counters ( values . data (), events . size ());
if ( retval != PAPI_OK ) {
print_papi_error ( retval );
}
// 測定結果の表示
auto & number_of_executed_instructions = values [ 0 ];
auto & number_of_floating_point_operation = values [ 1 ];
auto elapsed_cpu_time_usec = cpu_time_end_usec - cpu_time_begin_usec ;
auto GFLOPS = 1.e6 * number_of_floating_point_operation / elapsed_cpu_time_usec ;
std :: cout << "number_of_floating_point_operation = " << number_of_floating_point_operation << " \n " ;
std :: cout << "elapsed_cpu_time_usec = " << elapsed_cpu_time_usec << " \n " ;
std :: cout << "GFLOPS = " << GFLOPS * 1.e-9 << " \n " ;
std :: cout << "number_of_executed_instructions = " << number_of_executed_instructions << " \n " ;
// カウンタをクリア
PAPI_stop_counters ( values . data (), events . size ());
}
number_of_floating_point_operation = 100000000
elapsed_cpu_time_usec = 264001
GFLOPS = 0.378786
number_of_executed_instructions = 1200002473プログラムの出力のうち、 number_of_floating_point_operation だけは何回測っても変わらないという素晴らしい性質を持つ。他の値は測るたびに変わってしまう。
参考
PAPI
