ダブルダブル：FPUを離れることなく31桁の精度を実現

もしdouble形式の15桁の精度以上の精度が必要になった場合、2つのdoubleを結合して1つの数値として扱うという巧妙なトリックがあります。これにより、通常のdoubleの約9倍のコストで、約31桁の10進数精度が得られます。ヒープ割り当てや依存関係がなく、この手法はパフォーマンス面でdoubleと任意精度ライブラリのほぼ中間地点に位置します。この記事では、この手法を可能にするエラーフリー変換を説明し、MPFRと比較してその性能を測定し、このトリックが限界に達する場所を示します。 誰も埋めないギャップ 浮動小数点型には多くの「フレーバー」があります。例えば、floatは約7桁、doubleは約15桁の10進数精度を持ち、どちらもほぼ無料で利用できます。任意精度ライブラリは必要なだけ多くの桁数を提供しますが、操作ごとに大きなコストがかかります。「十分ではない」と「桁違いに遅い」の間にはギャップがあり、私はマンデルブロ集合の深いズーム中にそのギャップに陥りました。このページの最上部にある画像はこのギャップを示しています。同じビューが2回レンダリングされており、左側はdoubleの精度が尽きたときのブロック状の画像、右側はダブルダブルでレンダリングされたシャープな画像です。 マンデルブロ集合は有名なフラクタルで、無限に繰り返されるらせんやミニコピーで満たされています。深いズームは文字通り精度を使い果たします。やがて、隣接する2つのピクセル位置が同じdoubleに「丸め」られ、画像はフラクタルの画像ではなくなり、数値形式の画像になり始めます。これは図1で見ることができます。 図1: double精度でレンダリングされた2つの深い場所：ピクセルのブロック全体が1つの座標を共有するため、同じ色になります。double間のギャップは数値の大きさとともに増加するため、実軸は左側（re -0.74, im 0.13, 幅広のブロック）で粗く、虚軸は右側（re -0.11, im 0.92, 背の高いブロック）で粗くなります。 「double以上の精度が必要だ」という場合の標準的な答えは、通常GMPやMPFRのような任意精度数学ライブラリを使用することです。必要な精度が無限である場合には、これが正しい答えです。しかし、それは簡単なステップではなく、すべての操作でコストを支払うことになります。 すべての値はヒープブロックです。値はデータ配列（リム）へのポインタであるため、それを生成するには割り当てが必要です。多くのラッパーが提供するような、各操作が新しい値を返すAPIは、すべての操作で割り当てを行います。 すべての操作はループです。加算または乗算、すべてがリム配列をウォークします：ロード、ストア、データ依存分岐。 .NETではさらに悪いです。GMPとMPFRはネイティブCライブラリなので、P/Invoke、プラットフォームごとのネイティブバイナリの配布、境界でのマーシャリングが必要になります。 そしてライセンスの問題があります：GMPとMPFRはLGPLであり、静的リンク、クローズドプラットフォーム、または会社のポリシーによっては、ベンチマーク結果に関わらず、ハードな拒否となる可能性があります。 私はこれらすべてをベンチマークしましたが、ギャップは驚くほど構造化されていることがわかりました。同等の精度で実行した場合、実際のカーネルではダブルダブルは通常のdoubleの約9倍のコストがかかり、同じ31桁を扱うMPFRはダブルダブルのさらに約9倍のコストがかかり、元のdoubleの約81倍になります。MPFRが操作ごとに結果を割り当てる場合、ダブルダブルからMPFRへのギャップは9倍ではなく40倍に近くなります。 トリック全体は次のとおりです。 2つのdoubleを別々に保持しますが、それらを1つの数値として扱います。 約31桁の10進数精度、ヒープなし、依存関係なし、ネイティブバイナリの配布なし。 完全な型とベンチマークされたカーネルは2つの付録にあります。 ダブルダブルのアイデア 浮動小数点型は、数値の大きさに関わらず、同じ数の有効桁数を保持します。指数は小数点以下の位置を決定するだけです。これらの2つの数値を考えてみましょう。 (1)A = 111222333444 (2)B = 0.555666777888 どちらも12桁の有効数字を持ち、doubleに収まります。それらの正確な合計は111222333444.555666777888で24桁ですが、doubleは約15桁しか保持できません。A + Bを評価すると、収まらない桁は丸められ、111222333444.55566が残ります。 ここで重要なアイデアがあります：加算を評価しないとしたらどうなるでしょうか？ AとBを並べて保持し、ペアを1つの数値として扱うことに同意すれば、一方には先頭の桁が、もう一方には残りが格納されます。 これが完全な表現です。 (3)x = xhi + xlowith the invariant that xhi is exactly what xhi + xlo rounds to as a double. The high part carries the value, the low part carries the rest (the "error" of the high part after addition). Every operation in the next section is just keeping this invariant true. コードリスト1：ダブルダブル表現全体。この記事の他のすべては、これらの2つのフィールド上の算術演算です。 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 public readonly partial struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> { private readonly double m_hi; private readonly double m_lo; /// <summary> /// Leading component; equals the value rounded to nearest double. /// </summary> public double Hi => m_hi; /// <summary> /// Trailing component (the rounding error of <see cref="Hi"/>). /// </summary> public double Lo => m_lo; } 数学に入る前に、命名に関する注意点です。IEEE 754（doubleを定義する標準）の128ビット型はbinary128と呼ばれ、一般にquadruple-precisionとして知られています。これは連続した113ビットの仮数部と、約1e4932に達する15ビットの指数部を持ちます。 ダブルダブルも128ビットのストレージを使用しますが、単に2つの通常のdoubleを隣接させたものであり、合計約106ビットの仮数部と、変わらない指数範囲（依然として約1e308でオーバーフロー）をもたらします。 [1] 図2：各フォーマットの仮数部のビット数。 エラーフリー変換 すべては1つの洞察に基づいています。 double演算操作が丸めるとき、丸められた量は別のdoubleとして回復できる量です。 2つのdoubleを加算すると、ハードウェアは真の合計を計算し、それが丸められます。丸められた結果をs = round(a + b)とします。丸められた部分は正確に(a + b) - sであり、この式は自己矛盾しているように見えます：正確な合計が必要ですが、正確な合計は私たちが持っていない唯一のものです。それにもかかわらず、いくつかの通常のdouble操作で残りを正確に回復できます。アルゴリズムはKnuthのtwo-sumであり、コードリスト2はそのすべてです。6回の加算と減算で、分岐、ビット操作、より広い型はありません。 コードリスト2：Knuthのtwo-sum：任意の2つのdoubleに対して、sum + error == a + bが正確に成り立ちます。 1 2 3 4 5 6 7 8 private static (double Sum, double Error) twoSum(double a, double b) { double sum = a + b; double bKept = sum - a; // the part of b that made it into the sum double aKept = sum - bKept; // and the part of a double aLost = a - aKept; double bLost = b - bKept; return (sum, aLost + bLost); } 私たちの例の数値AとBで、このアルゴリズムの各行を見てみましょう。 コードリスト3：2つの例の数値に対するKnuthのtwo-sum。 1 2 3 4 5 6 7 8 9 10 11 12 13 A = 111222333444 B = 0.555666777888 A + B, exactly = 111222333444.555666777888 sum = A + B = 111222333444.5556640625 <- what the hardware hands back bKept = sum - A = 0.5556640625 aKept = sum - bKept = 111222333444 aLost = A - aKept = 0 bLost = B - bKept = 0.000002715388 error = aLost + bLost = 0.000002715388 sum + error = 111222333444.555666777888 <- exact, to the last digit doubleの合計は、正確な合計の先頭17桁を保持し、残りをドロップしました。これは予想通りです。[2] 注目すべき行はbLostです。これはドロップされたものの近似ではなく、ドロップされたものそのものであり、常にdoubleに収まります。なぜなら、カットオフ以下のすべてが小さい方のオペランドから来ているためであり、小さい方のオペランドは53ビットを超えることはなかったからです。極端な場合、オペランドが重なり合わないほど離れていると、合計は大きい方に丸められ、残りは小さい方全体になります。 これにより、最後の行がこの型全体が構築されているプロパティになります：(sum, error) は、最後の桁まで正確に A + B です。正確さは、配置によるものではなく、