RealLib ではじめる誤差ゼロ実数計算
RealLib のソースコード読みを始めるはずだったんですが、なんだか全然進んでないので適当なまとめエントリでお茶を濁します！ RealLib が普通にかっこよすぎるので紹介しまくりたくなりましたので紹介記事です。
実数計算と誤差
たいていのプログラミング言語の「実数 = 浮動小数点数」の計算には「誤差」があります。たとえばPythonのばあい：
Python 2.5 (r25:51908, Sep 19 2006, 09:52:17) [MSC v.1310 32 bit (Intel)] on win32 Type "help", "copyright", "credits" or "license" for more information. >>> 0.1 + 0.1 + 0.1 - 0.3 5.5511151231257827e-017
0.1 を 3 回足しても 0.3 にはなりません。少しだけずれちゃってます。
この例のばあい、誤差の原因は、「浮動小数点数」が2進数で表現されていることです。"0.1" は2進数で書くと無限小数なので、途中で打ち切って記憶されてしまい、正確な "0.1" とは少しずれた値になってしまうのです。
それでは困る場合のために、多くのプログラミング言語には、10進数で実数を記憶する「BigDecimal」のようなライブラリが用意されています。またまた Python のばあい：
>>> Decimal("0.1") + Decimal("0.1") + Decimal("0.1") - Decimal("0.3") Decimal("0.0")
ピッタリ "0.0" になりました！
では、Decimal のようなライブラリを使えば、実数計算の誤差は完全に消えてなくなるのでしょうか？ 答えは否。
>>> Decimal(8) / Decimal(7) * Decimal(7) Decimal("8.000000000000000000000000001")
"8/7" は10進数では無限小数になってしまうので、Decimalでも途中で打ち切られてしまいます。その結果、7倍しても元の8には戻らず、やっぱり微妙にずれます。Pythonのデフォルトでは途中の打ち切りは28ケタですが、これを調整してどれだけ長めにとっても、
>>> getcontext().prec = 80 >>> Decimal(8) / Decimal(7) * Decimal(7) Decimal("7.999999999999999999999999999999999999999999999999999999999999999999999 9999999997") >>> getcontext().prec = 300 >>> Decimal(8) / Decimal(7) * Decimal(7) Decimal("8.000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000002")
ずれるものはずれます。「7進数で表現するライブラリを導入すれば・・・！！」と思わず考えてしまいますが、何進数で表現しようとも、有理数で近似している限りは、たとえば sqrt や sin のような無理数を返す計算を入れたとたんに、ちょっとずれてしまうことは避けられません。（巧く実装しない限り）計算を繰り返すたびにこの誤差は少しずつ積もっていきます。
※ CPUの機能を使う浮動小数点演算と違って、BigDecimal 系のライブラリでは、打ち切る桁数を必要なだけいくらでも長く設定することができます。なので、もちろん、誤差が積もり積もっても問題ないくらい長めに桁数を取っておけば、この誤差は実用上は問題にならないで済みます。けど、どのくらいの誤差が発生しうるかをいちいち考えるのって大変です・・・よね？わたしは大変です。そんなこと気にしたくない！
Exact Real Arithmetic
・・・ということで、「誤差ゼロ実数計算」の登場です。ご紹介する RealLib は C++ のライブラリで、Real というクラスを提供しています。四則演算等々は定義されているので、普通の数値型と思って使うことができます。コード例：
#include <iostream> #include <iomanip> using namespace std; #include <Real.h> using namespace RealLib; int main() { InitializeRealLib(); { cout << setiosflags(ios::fixed); Real r = Real("8") / Real("7") * Real("7"); cout << setprecision(28) << r << endl; // 28桁目まで表示 cout << setprecision(200) << r << endl; // 200桁目まで表示 cout << setprecision(3000) << r << endl; // 3000桁目まで表示 } FinalizeRealLib(); }
実行結果：
> rltest 8.0000000000000000000000000000 8.000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000 8.000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 ...(中略)... 00000000000000000000000000000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000
たとえ何千桁まで表示しようとも、r はピッタリ 8 です。ずれません。ずれないのは有理数演算に限りません。
#include <iostream> #include <iomanip> using namespace std; #include <Real.h> using namespace RealLib; int main() { InitializeRealLib(); { cout << setiosflags(ios::fixed); Real s = cos( Pi/3 - Pi/2 ) * 2; // cos(-π/6)*2 = √3 Real ss = s * sqrt(Real("3")); // √3*√3 = 3 cout << setprecision(700) << s << endl << ss << endl; // 700桁表示 } FinalizeRealLib(); }
円周率やコサインや平方根の混ざった計算でも・・・
> rltest2 1.732050807568877293527446341505872366942805253810380628055806979451933016908800 03708114618675724857567562614141540670302996994509499895247881165551209437364852 80932319023055820679748201010846749232650153123432669033228866506722546689218379 71227047131660367861588019049986537379859389467650347506576050756618348129606100 94760218719032508314582952395983299778982450828871446383291734722416398458785539 76679580638183536661108431737808943783161020883055249016700235207111442886959909 56365797087168498072899493296484283020786408603988738697537582317317831395992983 00783870287705391336956331210370726401924910676823119928837564114142201674275210 23729942708310598984594759876642888977961478379583902288548529 3.000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000
ピッタリ 3 になります。√3 の方は表示を700桁で切った都合でピッタリ √3 ではありませんが（というか、ピッタリの√3を有限小数で表示するのは無理）、表示の量を増やせばいくらでも正確に√3に近い値が表示されます。Piやcosを使って遠回りに計算したからといって誤差が乗ることもないです。
使い方
適当に VC++ か gcc でビルドして適当に使うと使えると思います（適当ですいません）。公開されてるソースの Real.h は一カ所バグってて、確かどれか operator/= の実装が掛け算になってましたので直してから使った方が良いと思います（適当）。
よくある関数はひととおり提供されています。一覧表：
- コンストラクタ
- 文字列から構築 Real("12.3456")
- doubleから構築 Real(0.1) // これだと当然「正確な」0.1にはならないので注意
- 四則演算
- 加算 Real + Real
- 減算 Real - Real
- 乗算 Real * Real, Real * int, int * Real, Real * double, double * Real
- 除算 Real / Real, Real / int, int / Real, Real / double, double / Real
- 加算代入 Real += Real
- 減算代入 Real -= Real
- 乗算代入 Real *= Real, Real *= int, Real *= double
- 除算代入 Real /= Real, Real /= int, Real /= double
- 逆数 recip(Real)
- 符号反転 - Real
- 定数
- Pi
- Ln2
- 数学関数
スピード
精度に限りがないということはさぞかし計算が大変で重いだろう・・・とも思えますが、意外と、十分使えるスピードが出ます。例：
#include <iostream> #include <iomanip> #include <boost/progress.hpp> using namespace std; #include <Real.h> using namespace RealLib; int main() { InitializeRealLib(); { cout << setiosflags(ios::fixed); cout << setprecision(1000); { boost::progress_timer t(cerr); Real e(0.0); Real q(1.0); for(int i=0; i<100000; ++i, q/=i) e += q; cout << e << endl; } { boost::progress_timer t(cerr); double e(0.0); double q(1.0); for(int i=0; i<100000; ++i, q/=i) e += q; cout << e << endl; } } FinalizeRealLib(); }
10万回ループ。精度1000桁まで表示。速度面ではCPUの浮動小数点演算には敵うべくもないですが、手元のマシン（PenM 1.7GHz, VC8 /Ox）で1秒くらいで計算完了しました。
2.718281828459045235360287471352662497757247093699959574966967627724076630353547 ...(中略)... 041718986106873969655212671546889570350354 1.11 s 2.718281828459045500000000000000000000000000000000000000000000000000000000000000 ...(中略)... 000000000000000000000000000000000000000000 0.01 s
実際は500回も回せば1000桁の精度には達しているので、そこで打ち切ると0.05秒くらいです。色々な条件や適した使い方などなど違いますのであまり意味のある比較になってないのですが、一応、PythonのDecimalで精度1000桁の500回ループを回すと20秒弱かかります。
>>> from time import * >>> getcontext().prec = 1000 >>> def calc_e(): ... t = clock() ... e = Decimal(0) ... q = Decimal(1) ... for i in range(1,500): ... e += q ... q /= i ... print e ... print clock() - t, "sec" ... >>> calc_e() 2.718281828459045235360287471352662497757247093699959574966967627724076630353547 ...(中略)... 04171898610687396965521267154688957035044 19.4357855538 sec
比較演算落とし穴
一個だけ注意しないと問題なところがあります。比較演算子は以下のものが提供されています・・・
- Real != Real
- Real < Real
- Real > Real
・・・が、これで完全に等しい値どうしを比較すると、無限ループになって帰ってきません。無限精度の中で値が違う桁を見つけたら結果を返すという処理が行われるので、等しい値だとひたすら深い桁の値を調べに行ったっきり戻ってこなくなっちゃいます。現状では、「違う」と確実に分かっている値どうしの大小比較でないと使うのは危険かも。
妥協して、精度を指定して「この桁まで同じだったら同じと見なす大小判定」のようなものが用意されてると便利かもと思ったんですが、どうなのかな。
次回予告
おひさしぶりです！
次から「Exact Real Arithmetic （完全精度実数計算？定番の訳語ってなんだろう？）」のお勉強をはじめます。
のコードを読みつつ
を読みつつ、という形になる予定です。
も、もしかしたら参考にします。できれば見様見真似で自前の実装もしてみたいです。
いわゆる BigDecimal（Oracle Technology Network for Java Developers | Oracle Technology Network | Oracle、プログラミング言語 Ruby リファレンスマニュアル、http://www.python.jp/doc/nightly/lib/module-decimal.html、など）との違いは、BigDecimal は最初に設定した有効桁数の範囲でのみ正確に計算ができるのに対して、ExactReal は本当に無限の精度で計算ができること、と理解しています。どういう演算ならどのくらいの速度でどんな感じに計算できるのか〜という感覚をつかめたらなあと思います。
レジスタマシンとスタックマシン
こちらで見かけた話題・・・
Lua の VM は 4.0 でのスタックマシンから 5.0 でレジスタマシンに変更されたらしく、その時の理由が Lua: papers の "The implementation of Lua 5.0" の 第 7 節 にいくつか書いてありました。そこで書いてあることには・・・、の前に Lua 5.0 の「レジスタ」の仕様を簡単にまとめるておくと、こんな風になってます。
富豪的です。その分、レジスタ割り当てを考える必要がないので、コンパイラのコード生成が大変／実装が難しい、ということはないです。なにせ構文木をそのまま深さ優先探索でたどるだけでコードが出てきます。引数の問題、状態の保存の処理は call と return 命令に全部集約されているので、直接スタックを操作する命令も不要です。とゆーか、ジャンプ時のスタック長の整合性とか考えなくていい分スタックマシンより実装楽かもという気もします。
話を戻して、"The implementation of Lua 5.0" に書いてあったレジスタマシン化の利点ですが、
レジスタマシンはコピーが少ない
"a = b + c" という式を考えると、'a','b','c' をそれぞれの変数のローカル変数番号としたとき、スタックマシンだと
getlocal 'b' # 'b'の値をスタックにpush getlocal 'c' # 'c'の値をスタックにpush add # スタックポインタ2つ減らして和を push setlocal 'a' # スタックトップの値をpopして'a'に書き込み
こう、都合 3 回か 4 回かのコピーが発生しますが、レジスタなら
add 'a' 'b' 'c'
1 回。そして、特に Lua の場合、変数の「値」は型を表すタグ(4Byte)とポインタかdoubleのunion(8Byte)で計12Byteもある構造体なので、あんまり何回もコピーしたくないという事情があるそうな。「ポインタの下位ビットを利用してタグを云々〜」みたいなトリッキーな技を避けたいと思うと動的な言語の VM は同じ問題を抱えがち。
レジスタマシンの命令は word align しやすい
要は命令を固定長にしやすい。スタックマシンだと、「コードサイズが小さい」利点を保つために命令が可変長になってしまいがち（add 命令は 1 Byte で済むけど、getlocal や jmp 命令、定数ロード命令等々は、ローカル変数インデックスやジャンプ距離のように余計な情報がいるので 1 Byte では足りない）
LuLu 0.05 Released!
http://lulu.luaforge.net/index.ja.html
リリースしました！よろしくお願いします！
もっとソースを綺麗に整頓してからluaforgeに上げようと思っていたんですが・・・なんだか飽きちゃって(^^;)、すっごい中途半端です。ごめんなさい。でも、このままずるずると何もしないよりはババーンと出しちゃった方がいいかなあ、ということで。このリリースをもってLuaソースコード勉強会のまとめということにします。ご愛読ありがとうございました（？）
LuLu の実装面で残っている課題は
- とりあえずソースを綺麗にする（最低、coroutine.yield のグローバル変数はなくす）
- メタテーブルをちゃんと実装
- エラー処理をちゃんと実装してエラーメッセージを読みやすくする
- 多値を扱う組み込み関数で多値の末尾がnilなときの動きが変なのを直す
- dofile や load 系関数の実装 （Yueliang と合わせてみる？）
辺りです。Lua熱が高まったら開発再開します。協力してくださる方いらっしゃいましたら大歓迎です。それ以前に、なぜかLuaForgeのSSHやCVSにアクセスできないでいるのでやり方教えてくださる方は大大大歓迎です！
おまけ：LuLu (4)
コルーチンと getfenv/setfenv 実装しました。あと、今回はインチキして loadfile 関数でホストのLuaにコンパイル処理をやってもらうことにしたので、luac の出力ではなく、直接 .lua のソースも実行できるようになっています。
lua lulu.lua /usr/local/lua/test/sieve.lua
みたいな感じで、Luaについてくるテストプログラム動かせます。かっこいい！
今残っている Lua 5.1 と比べて足りてない部分、変な部分は
- 未実装：メタテーブル
- 未実装：debugモジュール
- 未実装：module/require/package
- 手抜き実装：pcall/cpcall
- 手抜き実装：dofile/load/loadfile/loadstring/string.dump
- 手抜き実装：string.dump
- 手抜き実装：function と thread が table として使えちゃうバグ（内部実装が漏れ出している）。tostringしたときにもtable 扱いされる
となっています。他のライブラリ（io, string, math）は、ホストLuaのライブラリに直接投げるという手抜き処理ですが、それによって特に問題は起きないはず。pcallやdofileは、それだとあんまり上手くいかないのでどうにかしないといけません。
次回は、この辺りなんとかしたりYueliangとくっつけたりするかもしないかも。The LuaJIT Project 勉強会に移行するかもしないかも。ではまたー。
環境
ここまでグローバル変数の扱いすっごく適当に見てたんですけど、
Foo = {x = 100} (function () setfenv(1, Foo) y = 2 z = x + y end)() print(Foo.x) -- 100 print(Foo.y) -- 2 print(Foo.z) -- 102
setfenv/getfenvというので、グローバル変数の参照するテーブルを差し替えられるんですね。モジュールを実装するのによく使われるみたいです。
さりげなくスルーしてしまったGETGLOBAL, SETGLOBAL命令の実装を見直してみると、こんな感じです。cl->envというテーブルに、値を読み書きしています。
case OP_GETGLOBAL: { TValue g; TValue *rb = KBx(i); sethvalue(L, &g, cl->env); lua_assert(ttisstring(rb)); Protect(luaV_gettable(L, &g, rb, ra)); continue; } case OP_SETGLOBAL: { TValue g; sethvalue(L, &g, cl->env); lua_assert(ttisstring(KBx(i))); Protect(luaV_settable(L, &g, KBx(i), ra)); continue; }
というわけで、getfenv/setfenvはcl->envを操作しているはず。cl->は現在実行中の関数オブジェクトなので、つまり関数ごとに環境を設定できるわけですね。
// lbaselib.c static int luaB_getfenv (lua_State *L) { getfunc(L, 1); if (lua_iscfunction(L, -1)) /* is a C function? */ lua_pushvalue(L, LUA_GLOBALSINDEX); /* return the thread's global env. */ else lua_getfenv(L, -1); return 1; } // lbaselib.c static void getfunc (lua_State *L, int opt) { if (lua_isfunction(L, 1)) lua_pushvalue(L, 1); else { lua_Debug ar; int level = opt ? luaL_optint(L, 1, 1) : luaL_checkint(L, 1); luaL_argcheck(L, level >= 0, 1, "level must be non-negative"); if (lua_getstack(L, level, &ar) == 0) luaL_argerror(L, 1, "invalid level"); lua_getinfo(L, "f", &ar); if (lua_isnil(L, -1)) luaL_error(L, "no function environment for tail call at level %d", level); } }
getfenvでは、関数オブジェクトそのものを指定するか、コールスタックの深さを整数で指定します。関数オブジェクトが直接来たときはいいとして、整数だった場合は、lua_getstackでコールスタックをたどって関数を取得しています。
LUA_API void lua_getfenv (lua_State *L, int idx) { ... switch (ttype(o)) { case LUA_TFUNCTION: sethvalue(L, L->top, clvalue(o)->c.env); break; ... case LUA_TTHREAD: setobj2s(L, L->top, gt(thvalue(o))); break; default: setnilvalue(L->top); break; } ... }
取得した関数オブジェクトのenvを返す感じで。あ、コルーチンを指定してコルーチンのグローバル環境を変えることも出来るんですね。あれ、そうでもないかも。調べ中。
setfenvも基本的に同じで、getfuncで対象の関数オブジェクトをとってきて、lua_setfenv
static int luaB_setfenv (lua_State *L) { luaL_checktype(L, 2, LUA_TTABLE); getfunc(L, 0); lua_pushvalue(L, 2); if (lua_isnumber(L, 1) && lua_tonumber(L, 1) == 0) { /* change environment of current thread */ lua_pushthread(L); lua_insert(L, -2); lua_setfenv(L, -2); return 0; } ... }
lua_setfenvはこう、envにテーブルを設定。
LUA_API int lua_setfenv (lua_State *L, int idx) { ... case LUA_TFUNCTION: clvalue(o)->c.env = hvalue(L->top - 1); break;
というわけでした。基本ライブラリのmodule関数もVMの言語実装部分的な仕組みとしてはgetfenvとsetfenvと同じように実装できそうです。
