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
• 数学関数
  • abs
  • sqrt
  • rsqrt (1 / sqrt)
  • log
  • exp
  • cos, sin, tan
  • acos, asin, atan, atan2
  • cosh, sinh, tanh

スピード

精度に限りがないということはさぞかし計算が大変で重いだろう・・・とも思えますが、意外と、十分使えるスピードが出ます。例：

#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

・・・が、これで完全に等しい値どうしを比較すると、無限ループになって帰ってきません。無限精度の中で値が違う桁を見つけたら結果を返すという処理が行われるので、等しい値だとひたすら深い桁の値を調べに行ったっきり戻ってこなくなっちゃいます。現状では、「違う」と確実に分かっている値どうしの大小比較でないと使うのは危険かも。

妥協して、精度を指定して「この桁まで同じだったら同じと見なす大小判定」のようなものが用意されてると便利かもと思ったんですが、どうなのかな。