RealLib ソースコード勉強会 (1:終)
RealLibの実装がどういう原理になってるのか解明しようシリーズです。シリーズといいながら、読み始めてみたらかなりシンプルだったので1回で終わらせちゃいます！
実数オブジェクトの構造
// Real.h class Real { private: RealObject *m_pObject; ... };
普段ユーザーが直接さわる Real クラスは、RealObject クラスを単純にラップしたものになっています。RealObject は、計算式の形をそのままツリー状に表現するためのクラスです。
// RealObject.h class RealObject { ... protected: virtual Encapsulation Evaluate() = 0; public: Encapsulation GetEstimate(); ... virtual RealObject* GetSibling(int index); ... };
GetSibling メソッドでツリーの子要素を取得するみたいです。Evaluate() を呼び出すと、「現在必要とされている十分な精度」を持った値を表現する Encapsulation というオブジェクトを返します。GetEstimate() は Evaluate() をラップして、無駄な再計算を減らすためのキャッシュを管理しているようです。
つまり
Real a("0.1"); Real b("0.3"); Real c = a + a + a - b;
このようなコードを書いた時点では、実際には足し算や引き算は行われないで
|
RealBinary(Minus)
| |
RealBinary(Plus) |
| | |
RealBinary(Plus) | |
| | | |
RealFromString("0.1") RealFromString("0.3")かわりに、こんなツリー構造（共通部分は共有するので、正確にはDAG構造）がcに記憶されます。実際に値が必要になった時点（表示するときなど）に再帰的にEvaluate()メソッドが呼ばれて、必要な精度の値が計算されます。ネタを明かすとシンプルですね。特に途中で正規化するようなこともないようです。そのままの式構造を持っています。
ツリーのノードにはこんな種類があります。
// RealObject.h // double型の値をそのまま返す class RealFromDouble : public RealObject { double m_Value; ... }; // 文字列表記された値をそのまま返す class RealFromString : public RealObject { char *m_pString; ... }; // 精度を指定するとその精度までの正確な文字列表記を返すオラクル（神託）関数 typedef const char* (*OracleFunction)(unsigned precision); class RealFromOracle : public RealObject { OracleFunction m_pOracle; ... }; // ０引数関数（≒ 定数、円周率πなど）を表す。精度を指定されると // その精度までの誤差入り値を返す class RealNullary : public RealObject { typedef Encapsulation (*FuncNullary)(unsigned prec, UserInt otherData); FuncNullary m_pFunc; ... }; // １引数関数（単項 - など）を表す。引数から計算できる限りの精度で計算する class RealUnary : public RealObject { typedef Encapsulation (*FuncUnary)(const Encapsulation &arg, UserInt otherData); FuncUnary m_pFunc; RealObject *m_pArg; ... }; // 2引数関数（+,-,*,/ など）を表す。引数から計算できる限りの精度で計算する class RealBinary : public RealObject { typedef Encapsulation (*FuncBinary)(const Encapsulation &left, const Encapsulation &right, UserInt otherData); FuncBinary m_pFunc; RealObject *m_pLeft; RealObject *m_pRight; };
例として、Real の足し算を行うと
// Real.cpp Real operator + (const Real &lhs, const Real &rhs) { return Real(Plus, lhs, rhs, 0); } Real::Real(Real::FuncBinary binfunc, const Real &lhs, const Real &rhs) { m_pObject = new RealBinary(binfunc, lhs.m_pObject, rhs.m_pObject, 0); ... }
RealBinary(Plus) ノードが作られて、これを Evaluate() すると
// Real.cpp Encapsulation RealBinary::Evaluate() { return m_pFunc(m_pLeft->GetEstimate(), m_pRight->GetEstimate(), m_iUserData); } // RealEncapsulation.h static inline Encapsulation Plus (const Encapsulation &lhs, const Encapsulation &rhs, UserInt i) { if (UsingMachinePrecision) return lhs.m_mach + rhs.m_mach; else return lhs.m_est + rhs.m_est; } // RealEstimate.cpp Estimate operator + (const Estimate &lhs, const Estimate &rhs) { LongFloat s(lhs.m_Value + rhs.m_Value); return Estimate(s, lhs.m_Error + rhs.m_Error + RoundingError(s, s.AdditionRoundingError())); }
誤差を計算に入れた足し算が実行されるようになります。
誤差入り値の表現
さきほどから何回かでてきた Encapsulation クラスは、MachineEstimate と Estimate というクラスの単純なラッパです。
// RealEncapsulation.h class Encapsulation { public: MachineEstimate m_mach; Estimate m_est;
必要な精度がマシンのdouble演算の範囲に収まる場合に高速化のために使われるのが MachineEstimate (double low, high のペア)で、それ以外の一般ケースに対応するのが Estimate です。使う側がこの差を気にしないで済むように差違を吸収するのが Encapsulation クラスらしいんですが、もっと良いクラス名があるんじゃないでしょうか・・・
Estimate クラスは、おおよその値を表す m_Value フィールドと、その値と本当の値の間の誤差見積もりである m_Error フィールドでできあがっています。
// RealEstimate.h class Estimate { private: LongFloat m_Value; ErrorEstimate m_Error; // error, |m_Value - real| < m_Error ...
LongFloat は、いわゆる普通の BigDecimal で、「指定された精度の範囲までは正確に実数を表現する」クラスです。特に変わった実装でもなかったので、詳細は省略。
// ErrorEstimate.h class ErrorEstimate { public: u32 m_Man; // >= 2^31 i32 m_Exp; // error = m_Man * 2 ^ (m_Exp - 32) ... };
誤差の表現である ErrorEstimate は仮数部と指数部による表現。この形式ということは、RealLib でも 1/2^(2^31) より小さい値は正確に表現しきれないようですね。
必要な精度の出し方
「計算式のツリー構造を作って、必要な時に必要な精度をもったEncapsulationをEvaluate()で計算する」という構成になっていることはわかりました。では、「必要な時」っていつでしょう。これには、３パターンあります。
- 値を表示するとき
- double へ変換するとき
- 大小比較するとき
他の計算はすべて、計算式のツリー構造をひたすら積み上げていくだけです。値の表示は、ようは桁数を指定しての文字列への変換なので、これとdoubleへの変換の時は、指定された必要精度になるようにEvaluate()します。
// Real.cpp double Real::AsDouble() const { Encapsulation::Computation comp; Encapsulation r(MakeGoodEstimate(m_pObject, 55, 1024)); // 55bit（0付近では1024bit）精度で計算 return r.weak_AsDouble(); // ErrorEstimateを無視してValueをdoubleに変換するメソッド } std::ostream& operator << (std::ostream& out, const Real &re) { Encapsulation::Computation comp; u32 pr = out.precision(); // 表示したい精度を取得 if (out.flags() & std::ios::scientific) ++pr; u32 bitswanted = u32(ceil(LOG_2_10 * (pr))); u32 pzbits = bitswanted; if (!(out.flags() & std::ios::fixed)) pzbits *= 2; Encapsulation r(MakeGoodEstimate(re.m_pObject, bitswanted, pzbits)); // その分計算 ...(以下略)... }
MakeGoodEstimate メソッドが肝心の処理を実行してます。
// Real.cpp Encapsulation MakeGoodEstimate(RealObject *pObj, i32 precision, i32 probzeroprec) { if (g_WorkingPrecision*32 <= precision + 32) InternalReset(precision/32 + 3);
g_WorkingPrecision が現在の計算精度を表すグローバル変数です。InternalResetは、細かいことを抜きにすると、g_WorkingPrecision を設定するメソッドです。このように、必要精度はグローバルにはグローバル変数で引き回されてる実装です。
u32 prec = g_WorkingPrecision; do { try { Encapsulation r(PerformEvaluation(pObj)); if (r.GetRelativeError() < precision && !IsProbableZero(r, probzeroprec)) { throw PrecisionException("MakeGoodEstimate()"); } return r;
PerformEvaluation は、式ツリーを再帰的にたどって、g_WorkingPrecision 精度で下から順に計算していくメソッドです。スタックオーバーフローしないように実際の再帰はスタック操作で頑張って書かれていますが、やることはそれだけなので詳細略。式の末端の精度を g_WorkingPrecision から始めるので誤差が積もって本来必要な精度に足りなくなることもありますが、足りればその値を返します。足りなければ例外を投げて・・・
} catch (PrecisionException e) { if (prec >= g_precMax) throw; else InternalReset(prec = min(g_precMax, prec * PRECMULTIPLIER / PRECDIVISOR)); } } while (1); }
・・・キャッチ。g_WorkingPrecisionを何倍かして、再計算。ループ。精度を上げる倍率は、現バージョンのソースコードでは2倍となっていました。これで末端精度を上げていけばそのうち最上位の精度も足りるようになります。なるほど。低精度の桁を何回も計算し直すことになるので遅そうですけど、幅優先と反復深化の差みたいなもので、実際にはこっちの方が効率的だったりするのかな。
大小比較は、RealLibでは「引き算して正か負かを判定」という実装になっています。つまるところ必要なのは、確実に正か負かわかる（＝誤差範囲にゼロが含まれなくなる）まで精度を上げる、という処理です。
// Real.h static inline bool operator > (const Real &lhs, const Real &rhs) { return (lhs - rhs).IsPositive(); } // Real.cpp bool Real::IsPositive() const { Encapsulation::Computation comp; Encapsulation r(MakeNonZero(m_pObject)); return r.IsPositive(); }
この処理は、MakeNonZero関数の担当。
Encapsulation MakeNonZero(RealObject *pObj)
{
u32 prec = g_WorkingPrecision;
do {
try {
Encapsulation r(PerformEvaluation(pObj));
if (!r.IsNonZero()) {
throw PrecisionException("MakeNonZero()");
}
return r;
} catch (PrecisionException e) {
if (prec >= g_precMax) throw;
else InternalReset(prec =
min(g_precMax, prec * PRECMULTIPLIER / PRECDIVISOR));
}
} while (1);
}
構造はMakeGoodEstimateとまったく同じですね。精度が必要な分に達したかの判定の代わりに、.IsNonZero() でゼロではなくなったかどうかを判定しています。あとは2倍ずつ精度をあげていくのみ。
個々の誤差入り演算
足し引き掛け割りは、概算値と最大誤差の間の計算によって見積もることができます。たとえば足し算は上で引用したように、誤作動しも足し算、掛け算はこんな感じ。
// RealEstimate.cpp Estimate operator * (const Estimate &lhs, const Estimate &rhs) { LongFloat r(lhs.m_Value * rhs.m_Value); ErrorEstimate e(lhs.m_Error * rhs.m_Error + lhs.m_Error * ErrorEstimate(rhs.m_Value) + rhs.m_Error * ErrorEstimate(lhs.m_Value)); return Estimate(r, e + RoundingError(r, r.MultiplicationRoundingError())); }
三角関数や指数関数、πなどの値は、テイラー展開やニュートン法などを使って、必要な精度に収束するまで近似を繰り返します。例えば π だと
// RealFuncs.cpp template <> Estimate pi(unsigned int prec) { ... prec = g_WorkingPrecision; Estimate v = recip(rpi(prec)); ... }
g_WorkingPrecision で必要な精度を引っ張ってきて、その精度で 1/π を計算して（πを直接算出するよりも楽）、逆数を求めて返す、という感じになってます。1/πの実装はBorwein iteration (8 - 13)という収束列だそうな。この級数は1ステップで精度が4倍に上がることが知られているらしく、その結果を使って必要な回数でループを止めます。
// RealFuncs.cpp Estimate rpi(unsigned int prec) ... for (u32 i = 32; i < prec*128; i = i*4) { ... } }
他の無限級数で計算される関数も止め方はだいたい同じようでした。
おしまい
他の同じ分野のライブラリも調べて回っているところなのですが、「末端精度を上げていけばそのうち最上位の精度も足りる」— ようするに、連続である、というのがこの手の実数計算の基本的な原理のひとつなみたいですね。（なので、ふつうの浮動小数点演算にはある floor/ceil/round のような不連続関数は提供されていないし、比較演算はちょうど等しいときにうまく動かなかったりする。）
RealLib の実装は、値を計算式の形で覚えておいて、精度が必要になったら、精度を上げて下から再計算、というものでした。逆向きに、値を遅延評価のリストとして計算しておいて、上の方で精度が必要になると遅延された評価が実行されて、下の方のリストの遅延評価をさらに進める・・・といった実装もあるようです。究極的にはほぼ同じことが起きる形になっていると思いますが、実装の雰囲気はかなり違ってきそうなので、そういうののソースも機会があれば読んでみたいところです。
