ひとり勉強会

ひとり楽しく勉強会

forall と -> (1.2)

いくつか細かい捕捉。

intro
introsと同じひっぺがし作業を1回だけやる
単数形にすると「1回だけやる」意味になります。introsだとやれるまでやってしまうので、今回の例だと4個全部はがしちゃいますが

Theorem Modus_ponens:
 forall (P: Prop),
   forall (Q: Prop),
     (P -> Q) -> P -> Q.
Proof.
intro.
intro.
intro.

3回だけで止めると

P: Prop
Q: Prop
H: P -> Q
------------------------------------------------
P -> Q

こうなるので、ここで exact H. という証明もありです。apply とか要らない分簡単かも。

intro X
X という名前をつけつつintroでひっぺがす
intros X1 X2 ... Xn
X1 X2 ... という名前をつけつつn個introsでひっぺがす
単にintroとかintrosとか打つと、Coqが適当に"H"やら"H0"と名前をつけてくれますが、自分で命名もできます。たぶん可読性のためにはこっちの方が良いです。

Theorem Modus_ponens:
 forall (P: Prop), forall (Q: Prop), (P -> Q) -> P -> Q.
Proof.
intros P Q P_naraba_Q.
exact P_naraba_Q.
Qed.

他にも、intro/intros にはいくつか賢いオプション指定ができるので、詳しくはマニュアル参照のこと。

assumption
仮定にあるどれかがゴールそのまんまだからそれでお願い!という証明

Theorem Modus_ponens:
 forall (P: Prop), forall (Q: Prop), (P -> Q) -> P -> Q.
Proof.
 intros.
 apply H.
 assumption.  (* exact H0. の代わり *)
Qed.

実は、"exact H" や "exact H0" のように名前を指定しなくても、そのくらい自動で探してこいよ!というとCoqさんちゃんと見つけてくれます。賢いですね。

賢いと言えば

auto
Coqさんならできる!頑張れよ!もっと熱くなれよ!!

Theorem Modus_ponens:
 forall (P: Prop), forall (Q: Prop), (P -> Q) -> P -> Q.
Proof.
 auto.
Qed.

実はCoqさん十分賢いのでこの程度の証明は、auto. って書くと勝手に全部やってくれるんですが、それだとあんまりCoq入門にならないので、この記事ではできるだけautoを使わないようにします。

forall と -> (1.1)

Modus Ponens

一番最初の例として

を証明してみよう。古代ギリシャから伝わる "Modus Ponens" とよばれる定理「"PならばQ" であり、しかも "P" ならば、それすなわち "Q" である」です。この定理を Coq で書くと、こうなります。

Theorem Modus_ponens:
  forall (P: Prop), forall (Q: Prop), (P -> Q) -> P -> Q.
  • 「forall ほげ ふが」は「どんなほげについても、ふがが成り立つ」と読 む。
  • 「ほげ -> ふが」は「ほげ ならば ふが」と読む。-> は右結合。
  • Anarchy Proof には forall (P Q: Prop), (P -> Q) -> P -> Q. と短く書いてありますが、略記法です。連続するforallはまとめて書けます。意味は同じです。

合わせると「どんな命題*1 P と Q を持ってきても、"P->Q" ならば、しかもさらに "P" ならば、"Q" である」。このThorem文の最後もピリオドで終わってますが、基本的に、Coqの文はすべてピリオド終端です。

証明

これを CoqIDE に打ち込んで、Ctrl+Alt+↓ を押すと、こうなる。Ctrl+Alt+↓ は、CoqIDE 上で1つ Coq のコマンドを実行するショートカットキーです。ちなみに Ctrl+Alt+↑ で1つ戻れます。

すると、右上に

1 subgoal
------------------------------------------------
forall (P Q: Prop), (P -> Q) -> P -> Q.

と出る。Coq さんが「対話証明モード」に入りました。「Theorem などと宣言したからには、この forall (P Q: Prop), (P -> Q) -> P -> Q. とやら、証明する覚悟はできているのであろうな????」という問いかけです。証明しましょう。

Theorem Modus_ponens:
  forall (P: Prop) (Q: Prop), (P -> Q) -> P -> Q.
Proof.

Proof. と打って Ctrl+Alt+↓。別に何も起こらないので書かなくていいんですけど、なんとなく気分が出るので Proof. と打ってから証明始めます。気分は重要です。

Theorem Modus_ponens:
 forall (P: Prop), forall (Q: Prop), (P -> Q) -> P -> Q.
Proof.
intros.

Ctrl+Alt+↓

対話証明モードに入ったら、「tactic」と呼ばれる、証明コマンドを入力することで証明をすすめます。"intros"は最重要tacticのひとつ。

intros
「"A -> B" を証明しろ」という状態を、「"A" が成り立つものと仮定して "B" を証明しろ」に変える。
"->" が沢山あっても一気に剥がす。「"A -> (B->C)" を証明しろ」を、「"A"と"B"が成り立つものとして、"C"を証明しろ」に。

まあ、"ならば" の定義みたいなものでしょうか。"A ならば B" を証明するには、「"A"が成り立つと仮定したら"B"が証明できる」ことを言えばよい。そういう方針で証明しますんでよろしくお願いします、とCoqさんに伝えるtacticが、introsです。ここら辺、言葉遊びをしているだけに聞こえるかもしれませんが、まあ先を見てみましょう。

intros. を実行すると、右上のCoqさんの要求がこのように変化します。

1 subgoal
P: Prop
Q: Prop
H: P -> Q
H0: P
------------------------------------------------
Q

横線の上が、いま仮定として使ってよいもの、下が証明すべきゴールです。Pという命題があって、「Qという命題があって、P->Qが成り立っていて、Pが成り立っている」 という仮定のもとで、Q を証明せよ、という状態になりました。あ、そうそう、言い忘れてましたが、intros さんは "->" だけでなくて "forall" も全部ひっぺがして仮定に置きます(というか、実は仕組み的には->もforallもCoqの中ではほとんど同じものだったり…)。

元々のゴールから

  1. "forall (P: Prop), "forall (Q: Prop), (P->Q) -> P -> Q" を証明せよ
  2. "P: Prop" があったと仮定して "forall (Q: Prop), (P->Q) -> P -> Q" を証明せよ
  3. "P: Prop" と "Q: Prop" があったと仮定して "(P->Q) -> P -> Q" を証明せよ
  4. "P: Prop" と "Q: Prop" があって "P->Q" が成り立つと仮定して "P -> Q" を証明せよ
  5. "P: Prop" と "Q: Prop" があって "P->Q" と "P" が成り立つと仮定して "Q" を証明せよ

まで、4個のひっぺがしが行われました。

さて、お次のtacticは…。

apply H.

Ctrl+Alt+↓。ゴールが変わります。

1 subgoal
P: Prop
Q: Prop
H: P -> Q
H0: P
------------------------------------------------
P

違いがわかりますでしょうか。さっきまで "Q" だったゴールが "P" になってます。

apply H
「"A->B"の形の仮定Hがあって、ゴールが"Bを証明せよ"」という状態を「"Aを証明せよ"」というゴールに変える
これも、別方向からみた "ならば" の定義みたいなもの。"A->B" を仮定してよいのならば、"A"さえ証明できれば、そこから仮定を使って(applyして)"B"はすぐに証明できます。なので「"A->B"の形の仮定Hがあって、ゴールが"Bを証明せよ"」という状態では、直接"B"じゃなくても"A"さえ証明できれば勝利なのです。そういう方針で証明しますんでよろしくお願いします、とCoqさんに伝えるtacticが、applyです。

exact H0.

Ctrl+Alt+↓。さっきの証明目標図をよくみると、仮定に"P"があって、ゴールも"P"です。「"P"が成り立つことを仮定して、"P"が成り立つことを証明せよ」と言われているわけです。おまえさん、それはもう、仮定そのまんま(exact)じゃん!!とCoqさんに伝えるtacticが、exactです。

はい、これで証明完了しました。ゴールが "Proof completed." になるはずです。
最後に、"Qed." と打ち込んで締めます。これは"Proof."と違って気分の問題ではなく、打ち込まないと証明が終わりません。

Qed.

Ctrl+Alt+↓。

証明完了!

Coqの入門記事を書く会 (1)

証明支援系Coqの遊び方の入門を書いてみるよ。すでに世の中に何個も入門記事ありますけど、増えて困ることはなかろう…ということで。まず、インストール方法については、id:yoshihiro503 による紹介記事

最近パソコンを買ったんだけど、使い方がよく解らない。という人のために、今日はCoqのインストールの方法を紹介しよう。

を参照のこと。現時点では最新の安定版の Coq 8.2 がおすすめ。

充足と恒偽のあいだ

あるいは、ランダムなSAT問題の答えを分ける不思議な数字、4.27。

調べたり実験したりしてたら面白くなってきたので、まとめておきます。

3-SAT?

3-SAT は、「「「変数か、変数を NOT した式」を3個 OR で繋げた式」を何個か AND で繋げた式」が与えられたときに、「式全体の値が true になるように変数の値を決めることができるか?」を判定する問題です。NP完全問題の代表例として有名です。たとえば

(x1 || x2 || x3) && (!x1 || !x2 || !x3) && (x1 || !x3 || !x4)

この論理式は、x1=true、x2=false、x3とx4は適当に決めれば全体は true になるので、答えは "SATISFIABLE (充足可能)" です。
一方で、ちょっと人工的ですが、

(x1 || x1 || x1) && (!x1 || !x1 || !x1)

これは、x1 を true にしても false にしても全体は false なので、"UNSATISFIABLE (充足不可能)" が答えになります。

ランダム生成

3-SAT問題の論理式をランダム生成した場合、答えが "SATISFIABLE" になる確率と、"UNSATISFIABLE" になる確率、どちらが高いでしょうか?

ちょっと考えると感覚的にはなんとなくわかると思うのですが、3-SAT の論理式の場合、条件を && で繋いでいくので、長い論理式であればあるほど制約がきつくなって "UNSATISFIABLE" になりやすくなります。

  • 式が短いほうが SATISFIABLE になりやすい
  • 式が長いほうが UNSATISFIABLE になりやすい

逆に、式が長くても、たとえば10000項の論理式であっても、その中に変数が20000種類含まれていたら、変数の値を調整できる部分が大きいので、巧いことやって式をtrueにする抜け道が見つけやすくなります。

  • 変数の数が少ない方が UNSATISFIABLE になりやすい
  • 変数の数が多い方が SATISFIABLE になりやすい

…と、傾向としては、こうなりそうです。

Phase Transition!

面白いのは、この傾向は、どうも「なりやすい」というレベルではなく、ある点を境にかなり急激に変化するらしい。この手の問題の "Phase Transition" と呼ばれているそうです。

  • (式の長さ / 変数の数) ≦ 4.0 だとほぼ確実に SATISFIABLE
  • (式の長さ / 変数の数) = 4.27 付近で半々
  • 4.5 ≦ (式の長さ / 変数の数) だとほぼ確実に UNSATISFIABLE

(※ここでは、"式の長さ" は && で繋ぐ単位である (? || ? || ?) の個数を指しています。)

実験してみました。変数の数20で式の長さ3.0〜6.0倍まで色々変えてみたもの(赤)と、変数の数50で式の長さを3.0〜6.0倍まで色々変えてみたもの(青)。解くのは PicoSAT に頑張ってもらいました。それぞれの点で100回ランダム生成した問題を解いてSATISFIABLEだった率を集計してます。

まだ問題サイズが小さすぎるのか、3.5〜5.0の辺りまであいまいな部分が広がっていますが、かなり傾向は見て取れます。"Experimental results on the crossover point in random 3-SAT" という論文に200変数の場合の実験結果が載っていて、かなり綺麗なグラフになっています。以下に引用します。

逆に言うと、4.27付近は「SATISFIABLEかUNSATISFIABLEか全く予想がつかない」ので、SATを解くアルゴリズムにとってとても難しいゾーンということになります。上に引用したグラフにありますが、ヒューリスティックを駆使して効率的にSATを解くアルゴリズムは、式長/変数数比がこのゾーンの外ならとても効率的にSATを解けるらしいです。この4.27ゾーンだけが鬼門。

実験だけではなく、"A new upper bound for 3-SAT" によると、問題サイズを大きくしていくと、理論的に3.52以下はほぼ確実にSATISFIABLE、4.4898以上はほぼ確実にUNSATISFIABLEとなることが証明されていて、この理論値の幅はじわじわと狭められていっているそうな。

SATソルバを作る会 (1)

適当にまったりと。今回は下準備みたいなものだけです。

入力

せっかくなのでデータ形式SAT Competitions の形式と揃えておこうかと思いました。出せるようなものはとても作れないですけど、比較する時にその方が便利そう。

まず、入力はConjunctive Normal Formの論理式で与えるとのこと。

(x1 or !x5 or x4) and (!x1 or x5 or x3 or x4) and (!x3 or !x4)

こういう、『「変数か変数の否定 (= リテラル)」をorで繋いだもの (= 選言節)』をandで繋いだもの、という形の式です。ファイル形式としては、DIMACS形式というので与えられるらしい。こんなテキストファイルです。

c
c コメント〜
c
p cnf 5 3
1 -5 4 0
-1 5 3 4 0
-3 -4 0
  • 最初に何行かコメント (c XXXXX)
  • 次にヘッダ的なもの (p cnf 変数の個数 選言節の数)
  • 以下、一行に一個の選言節。スペース区切りの整数列で、
    • 正の数はnotがついてない変数
    • 負の数はnotがついてる変数
    • ゼロは選言節の終わり

出力

出力は、論理式を true にすることができる場合は

s SATISFIABLE
v 1 -2 3 -4 -5 0

のようにして、たとえば x1,x3 をtrueにして x2,x4,x5 をfalseにすればOK、という例を標準出力に出します。複数のtrueにしかたがある場合、どれか1つ出せばよいです。同時に、実行ファイルの終了コードとして 10 を返します。

出力は、論理式を true にすることができない場合は、

s UNSATISFIABLE

だけでいい。実行ファイルの終了コードとして 20 を返す。

コンテストのように時間制限がある場合、時間切れ等でどちらともわからなかった場合は、

s UNKNOWN

で、終了コードとして 0 を返す。とりあえずこれは自分には関係ないかな。

Take 1 : 完全に全探索

まずはウォームアップとして、全ての変数割り当てパターンを総当たりで試してみるソルバを書いてみました。今のところD言語を使っています。

import std.stdio;
import std.stream, std.cstream;
import std.string;
import std.conv;

class SAT
{
private:
  alias int Literal; /// 非ゼロ整数。+n は xn、-n は !xn を表す
  int  varOf( Literal lit )      { return (lit>0 ? lit-1 : -lit-1); }
  bool isPositive( Literal lit ) { return (lit>0); }

  alias Literal[] Clause; /// [-1,+2,+3] が !x1|x2|x3 を表す
  alias Clause[]     CNF; /// [[-1,+2],[+3,-2]] で (!x1|x2)&(x3|!x2) を表す

private:
  int numVar; /// 変数の個数
  CNF cnf;    /// 論理式

public:
  /// DIMACS format の入力ストリームから適当に読み込み
  this( InputStream sin )
  {
    foreach(char[] s; sin)
      if( s.length && s[0]!='c' )
        if( s[0] == 'p' )
          numVar = to!(int)( (s.idup.split)[2] );
        else
          cnf ~= to!(int[])( (s.idup.split)[0..$-1] );
  }

  /// ビットベクトルで変数への値割り当てを与えて評価
  bool eval(BitVec)( BitVec bv )
  {
    bool be = true;
    foreach(clause; cnf)
    {
      bool bc = false;
      foreach(lit; clause)
      {
        static if( is(bv[0]) )
          bool bl = bv[varOf(lit)];
        else
          bool bl = !!((bv>>varOf(lit))&1);
        bc |= isPositive(lit) ? bl : !bl;
      }
      be &= bc;
    }
    return be;
  }
}

bool solve( SAT s )
{
  long MASK = (1L << s.numVar) - 1;

  B_LOOP:
  for(long B=0; B<=MASK; ++B)
    if( s.eval(B) )
    {
      // Satisfied
      writeln("s SATISFIABLE");
      write("v");
      for(int i=0; i<s.numVar; ++i)
        write(" ", (B>>i)&1 ? i+1 : -i-1);
      writeln(" ", 0);
      return true;
    }

  // Not satisfied
  writeln("s UNSATISFIABLE");
  return false;
}

int main()
{
  auto s = new SAT(din);
  return solve(s) ? 10 : 20;
}

さて、これが最低ラインのパフォーマンスと言うことになります。どのくらいのスピードかなー。

20変数85節の3-CNF(各選言説のリテラル数が全部3のもの)ランダムデータを作って試してみました。色々試してみたところ、節数をこれより減らすとSATISFIABLEだらけ、これより増やすとUNSATISFIABLEだらけになってしまうようなので。比較対象はPicoSATです。

<picosat>
s SATISFIABLE
v -1 2 3 -4 -5 6 7 8 9 -10 11 12 13 -14 15 -16 17 -18 19 20 0
[10] : 99 ms

<hzkrsat>
s SATISFIABLE
v 1 2 -3 -4 5 -6 7 8 -9 10 -11 12 -13 14 -15 16 -17 -18 19 -20 0
[10] : 5977 ms

といったところ。もちろん、もっと大きい入力を喰わせるとどんどん差は開いていきました。25変数105節だと…

<picosat>
s SATISFIABLE
v 1 -2 3 -4 5 6 7 -8 -9 -10 11 -12 -13 -14 -15 -16 -17 18 19 20 -21 22 -23 -24
v 25 0
[10] : 105 ms

<hzkrsat>
s SATISFIABLE
v -1 2 -3 4 5 6 -7 -8 -9 10 -11 -12 -13 -14 -15 -16 17 18 19 -20 21 -22 -23 24 -
25 0
[10] : 222315 ms

というわけで、次回から色々試して改善していけたらいいなーと思います。

雑談

そういえば、変数と節数の比が4.4898を超えるとほぼ確実にUNSATISFIABLEという論文を前に読んで「へー!」と思っていたのを思い出しました。確かに自分で手で実験してみるとその辺りに境がありますね。面白かったです。

ググってみると、イマドキのSATソルバは基本的に DPLL という枠組みで解いているらしい。次はこの辺りをちょっと調べてみようかと思います。予定は未定です。

次回予告など

すっごく遅いけどあけましておめでとうございます!

忙しくて他のことをする余裕がないときほど突然に他のことをやりたくなるものです。でも、いざ忙しさのピークが過ぎる頃には情熱が冷めていたりして結局何もやらないまま過ぎてしまう。

これは良くないので、やるぞと世界に向けて宣言しておいて自分を追い込むメソッドを発動しようと思います。

というわけで来週末から『新しい SAT Solver を作る会』やります。

SAT( 充足可能性問題 - Wikipedia )というのはNP完全問題‐効率的に解くのが難しいとされている問題の代表格ですが、実際には、いろいろなヒューリスティックスで挑むことで、計算量のオーダから想像するよりも遙かに速いスピードで解くことができるらしい。この前「いいニュースと悪いニュースがある。悪いニュースは、この問題はNP完全ということだ。いいニュースは、簡単にSATに還元できるから実際は効率的に解けるということだ」というセリフを聞いてすごぃ!と思ったので、影響されてます。

とりあえず自分で考えるだけ考えた実装したあとは、天下一SAT武道会、じゃなくて SAT Competitions の上位のソースなどを読みながら遊んでみたいと思ってます。

あと SICMグループ ができてるのをみて思わず参加してしまったので、こっちも再開したいな。。。

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;
    ...
};
// 0引数関数(≒ 定数、円周率πなど)を表す。精度を指定されると
// その精度までの誤差入り値を返す
class RealNullary : public RealObject {
    typedef Encapsulation (*FuncNullary)(unsigned prec, UserInt otherData);
    FuncNullary m_pFunc;
    ...
};
// 1引数関数(単項 - など)を表す。引数から計算できる限りの精度で計算する
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()で計算する」という構成になっていることはわかりました。では、「必要な時」っていつでしょう。これには、3パターンあります。

  • 値を表示するとき
  • 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 の実装は、値を計算式の形で覚えておいて、精度が必要になったら、精度を上げて下から再計算、というものでした。逆向きに、値を遅延評価のリストとして計算しておいて、上の方で精度が必要になると遅延された評価が実行されて、下の方のリストの遅延評価をさらに進める・・・といった実装もあるようです。究極的にはほぼ同じことが起きる形になっていると思いますが、実装の雰囲気はかなり違ってきそうなので、そういうののソースも機会があれば読んでみたいところです。