付録E: 高精度演算（double-double型とxprec-rs）

標準の64ビット浮動小数点数 f64 は、ほとんどの科学技術計算において十分な精度と広い表現範囲を提供します。しかし、浮動小数点演算と誤差で見たように、浮動小数点演算には常に丸め誤差や桁落ちといった問題が伴います。特に、下記のようなケースでは f64 の精度では不十分となり、計算結果に大きな誤差が生じることがあります。

条件数の悪い問題: 行列計算などで、入力の小さな変化が出力に大きな影響を与える問題。
カオス系のシミュレーション: 初期値鋭敏性を持つため、初期の微小な誤差が時間とともに指数関数的に増大する。
長時間のシミュレーション: 繰り返し計算の過程で丸め誤差が累積し、無視できない大きさになる。

このような場合、標準の f64 よりも高い精度を持つ数値型が必要になります。本節では、その一つの解決策である double-double演算 と、それをRustで実現するためのクレート xprec を紹介します。

double-double演算の仕組み

double-double演算は、一つの数値を二つの f64 の和として表現することで、実質的により高い精度を実現する手法です。

表現: ある数値を、のように表現します。
- は数値の主要部分（上位ワード）を表す f64 です。
- はで表現しきれなかった微小な部分（下位ワード）を表す f64 です。

この方法により、f64 の仮数部が53ビットであるのに対し、double-double型では実質的に106ビットの仮数部を持つことになり、有効桁数を約2倍（10進数で約15桁から約31桁へ）に増やすことができます。

これは、IEEE 754で規格化されている四倍精度浮動小数点数とは異なり、既存の倍精度浮動小数点数の演算機能を組み合わせてソフトウェア的に高精度を実現する手法です。そのため、特別なハードウェアを必要としない利点があります。

`xprec`クレートによる高精度計算

xprec は、このdouble-double演算を効率的に行うためのRustクレートです。f64 とほぼ同じように扱える Df64 型を提供します。

`xprec`の導入

まず、Cargo.toml の [dependencies] セクションに xprec を追加します。

[dependencies]
xprec = "0.2"

`Df64`型の基本的な使い方

xprec の中心となるのは Df64 型です。f64 と同様の演算子（+, -, *, /）が実装されており、直感的に使用することができます。

use xprec::Df64;

fn main() {
    // f64からDf64への変換
    let a = Df64::from(1.0);
    let b = Df64::from(2.0);

    // 基本的な四則演算
    let sum = a + b;
    let diff = a - b;
    let prod = a * b;
    let quot = a / b;

    // Df64はDisplayトレイトを実装している
    println!("a + b = {}", sum);
    println!("a - b = {}", diff);
    println!("a * b = {}", prod);
    println!("a / b = {}", quot);

    // Df64からf64への変換（上位ワードを取得）
    let sum_f64: f64 = sum.hi();
    println!("Sum as f64: {}", sum_f64);
}

Note

Df64::from(x) の動作について

Df64::from(x) は f64 の値をそのまま上位ワードとして使用し、下位ワードを0に設定します。そのため、f64 の時点で既に丸め誤差が含まれている値（例: 0.1）は、Df64 に変換しても誤差は解消されません。Df64 の真価は、高精度な値同士の演算において誤差の累積を抑制できる点にあります。

実例: 総和計算における誤差累積の軽減

Df64 の有用性を示す具体例として、Basel問題（バーゼル問題）を取り上げます。これは、自然数の逆数の二乗和を求める問題で、解析解がであることが知られています。

この級数を有限項で打ち切って計算すると、f64 では項数が増えるにつれて丸め誤差が累積し、精度が悪化します。一方、Df64 を使用すると、誤差の累積を大幅に抑制できます。

use xprec::Df64;

fn main() {
    let n: u64 = 100_000_000; // 1億項
    let exact = std::f64::consts::PI.powi(2) / 6.0;

    let mut sum_f64: f64 = 0.0;
    let mut sum_df64: Df64 = Df64::ZERO;

    for i in 1..=n {
        let i_f64 = i as f64;
        let term = 1.0 / (i_f64 * i_f64);
        sum_f64 += term;
        sum_df64 = sum_df64 + Df64::from(term);
    }

    // 級数の打ち切り誤差を補正（1/N で近似）
    let correction = 1.0 / (n as f64);
    sum_f64 += correction;
    sum_df64 = sum_df64 + Df64::from(correction);

    println!("解析解:       {:.16}", exact);
    println!("f64  の結果:  {:.16}", sum_f64);
    println!("Df64 の結果:  {:.16}", sum_df64.hi());
    println!("f64  の誤差:  {:.2e}", (sum_f64 - exact).abs());
    println!("Df64 の誤差:  {:.2e}", (sum_df64.hi() - exact).abs());
}

実行結果（N = 10^8 の場合）:

解析解:       1.6449340668482264
f64  の結果:  1.6449340678345750
Df64 の結果:  1.6449340668482264
f64  の誤差:  9.86e-10
Df64 の誤差:  0.00e+00

この結果から、1億回の加算を行った場合、f64 では約10桁目に誤差が現れる（誤差 ≈ 10^-10）のに対し、Df64 では誤差がほぼ0に抑えられていることがわかります。これは、長時間のシミュレーションや大規模な数値計算において、Df64 が非常に有効であることを示しています。

数学関数の利用

xprec は simba クレートのトレイトを実装しており、sqrt, sin, cos, exp などの数学関数や、π などの定数を利用できます。これらを使用するには、対応するトレイトをインポートする必要があります。

use xprec::Df64;
use simba::scalar::{ComplexField, RealField};

fn main() {
    let x = Df64::from(2.0);

    // 平方根
    let sqrt_x = x.sqrt();
    println!("sqrt(2) = {}", sqrt_x);

    // べき乗
    let x_cubed = x.powi(3);
    println!("2^3 = {}", x_cubed);

    // 高精度のπ
    let pi: Df64 = Df64::pi();
    println!("π (Df64) = {}", pi);
    println!("π (f64)  = {}", std::f64::consts::PI);

    // 三角関数
    let sin_pi = pi.sin();
    println!("sin(π) = {}", sin_pi);
}

Important

トレイトのインポートについて

数学関数を使用するには simba クレートのトレイトをインポートする必要があります。Cargo.toml に simba を追加するか、xprec が依存している simba を利用してください。
[dependencies]
xprec = "0.2"
simba = "0.9"

利用上の注意点

double-double演算は強力ですが、万能ではありません。利用する際には以下の点に注意が必要です。

パフォーマンス: Df64 の演算は、複数の f64 演算を組み合わせて実行されるため、ネイティブの f64 演算よりも遅くなります（およそ15倍程度）。そのため、プログラム全体を Df64 に置き換えるのではなく、精度が特に要求される部分に限定して使用するのが効果的です。
任意精度ではない: double-doubleは、f64 の約2倍の精度を提供する固定精度の数値型です。さらに高い精度（例えば100桁など）が必要な場合は、rug や malachite といった任意精度演算ライブラリを検討する必要があります。

計算物理学において、ほとんどの場面では f64 で十分です。しかし、その限界と、xprec のような高精度計算の選択肢があることを知っておくことは、より信頼性の高いシミュレーションを行う上で非常に重要です。

Keyboard shortcuts

Rust計算物理学チュートリアル