浮動小数点演算と誤差

計算物理学では、連続的な物理量をコンピュータ上で離散的に扱います。その際、実数を表現するために浮動小数点数が用いられますが、この表現は有限のビット数で行われるため、本質的に近似であり、誤差の発生は避けられません。

本節では、この浮動小数点数に起因する誤差について解説し、それが数値計算の結果に与える影響と、その対策について学びます。誤差の存在を正しく理解し、適切に対処することは、信頼性の高いシミュレーションを行う上で極めて重要です。

Rustにおける浮動小数点数

Rustの標準的な浮動小数点数型には、f32とf64の2種類があります。 どちらも数学的な実数そのものではなく、有限のビット数で実数を近似した値です。

f32はf64よりメモリ使用量が少ないため、大きな配列ではメモリ使用量や メモリ帯域の面で有利になることがあります。ただし、常に高速とは限りません。 速度は計算内容、ハードウェア、メモリアクセスのパターンに依存します。

Important

本書ではf64を標準とします

計算物理学の多くのシミュレーションでは、計算過程で発生する微小な誤差が積み重なり、最終的な結果に大きな影響を与えることがあります。そのため、本書では特に断りがない限り、より高い精度を持つf64型を標準として使用します。

複素数

Fourier 変換や量子力学では、複素数を使います。 Rustの標準ライブラリには複素数型が含まれていないため、 本書では外部クレート num-complex を使います。

Cargo.toml に以下を追加します。

[dependencies]
num-complex = "0.4"

num-complex では、実部と虚部が f64 の複素数型として Complex64 を使えます。

use num_complex::Complex64;

fn main() {
    let z = Complex64::new(1.0, 2.0); // 1 + 2i
    let w = Complex64::new(0.5, -1.0);

    println!("z + w = {}", z + w);
    println!("|z| = {}", z.norm());
}

複素数でも、実部と虚部は浮動小数点数です。 そのため、丸め誤差や比較の注意点は f64 と同じように現れます。 複素数の詳しい使い方は、Fourier 解析や量子力学の章で扱います。

浮動小数点数に起因する誤差

浮動小数点数を用いる際に発生する代表的な誤差を3つ紹介します。

1. 丸め誤差 (Rounding Error)

コンピュータが2進法で数値を扱うことに起因する誤差です。10進法では有限桁で表現できる小数（例: 0.1）が、2進法では無限小数となり、有限のビット数で正確に表現できない場合があります。その結果、最も近い表現可能な値に「丸め」られることで生じるのが丸め誤差です。

この影響を確認するために、簡単なコードを実行してみましょう。

fn main() {
    let a: f64 = 0.1;
    let b: f64 = 0.2;
    let sum = a + b;
    let expected: f64 = 0.3;

    println!("0.1 + 0.2 = {}", sum);
    println!("期待値     = {}", expected);
    println!("両者は等しいか？ => {}", sum == expected);

    // 内部的な表現を確認
    println!("実際の値 (内部表現): {:.20}", sum);
}

このコードを実行すると、以下のような結果が得られます。

0.1 + 0.2 = 0.30000000000000004
期待値     = 0.3
両者は等しいか？ => false
実際の値 (内部表現): 0.30000000000000004441

0.1 + 0.2 の計算結果が、期待する 0.3 とは微妙に異なる値になっていることがわかります。これは、0.1 と 0.2 が2進数で正確に表現できず、それぞれが近似値として格納されているために発生します。

Warning

浮動小数点数の直接比較は避ける

丸め誤差のため、== 演算子による浮動小数点数の直接比較は予期せぬ結果を招く可能性があります。計算結果が特定の値と等しくなることを期待するようなロジックは避けるべきです。比較の方法については後述します。

2. 桁落ち (Cancellation Error)

桁落ちとは、値がほぼ等しい2つの数値の差を計算した際に、有効桁数が大幅に失われてしまう現象です。

例えば、二次方程式 の解の公式を考えます。

ここで、 の場合、 となります。もし であれば、一方の解 の分子において、 と がほぼ等しい値となり、桁落ちが発生します。

以下のコードは、 の解をナイーブな解の公式で計算する例です。

fn main() {
    let a: f64 = 1.0;
    let b: f64 = 1.0e8; // 1億 (大きな値にすることで桁落ちが顕著になる)
    let c: f64 = 1.0;

    let d = (b * b - 4.0 * a * c).sqrt();

    // -b と d は非常に近い値になる (絶対値が近い)
    println!("-b = {}", -b);
    println!(" d = {}", d);

    // 桁落ちが発生する計算
    let x1_naive = (-b + d) / (2.0 * a);

    // もう一方の解
    let x2 = (-b - d) / (2.0 * a);

    println!("\nナイーブな計算:");
    println!("x1 = {:.17}", x1_naive);
    println!("x2 = {:.17}", x2);

    // 桁落ちを回避する工夫
    // 解と係数の関係 x1 * x2 = c/a を利用する
    let x1_stable = c / (a * x2);

    println!("\n安定な計算:");
    println!("x1 = {:.17}", x1_stable);

    // 相対誤差を計算
    let relative_error = (x1_naive - x1_stable).abs() / x1_stable.abs();
    println!(
        "\n相対誤差: {:.2e} ({:.1}%)",
        relative_error,
        relative_error * 100.0
    );
}

実行結果：

-b = -100000000
 d = 99999999.99999999

ナイーブな計算:
x1 = -0.00000000745058060
x2 = -100000000.00000000000000000

安定な計算:
x1 = -0.00000001000000000

相対誤差: 2.55e-1 (25.5%)

ナイーブな計算で得られた x1 は、本来の値（約 ）と比較して精度が大幅に低下しています。一方、解と係数の関係 を用いて と計算することで、桁落ちを回避し、より精度の高い解を得ることができました。

Note

なぜ x2 は桁落ちしないのか

x2 = (-b - d) / (2a) において、 のとき は負、 は正であるため、 は「負の値からさらに正の値を引く」計算となり、絶対値が増加する方向の演算です。桁落ちは「ほぼ等しい値の減算」で発生するため、x2 の計算では問題が生じません。このように、減算を回避する数式の変形が桁落ち対策の基本となります。

このように、計算順序やアルゴリズムを工夫することで、桁落ちを回避できる場合があります。

3. 情報落ち (Loss of Significance)

絶対値が非常に大きい数値と非常に小さい数値を足し算する際に、小さい数値が結果に反映されず、情報が失われてしまう現象を情報落ちと呼びます。

fn main() {
    let large_number: f64 = 1.0e16;
    let small_number: f64 = 1.0;

    let result = large_number + small_number - large_number;

    println!("(1.0e16 + 1.0) - 1.0e16 = {}", result);
    println!("期待値 = 1.0");
}

このプログラムの実行結果は 0.0 となります。1.0e16 + 1.0 の計算において、1.0 は f64 が表現できる有効桁数の範囲外となり、無視されてしまったためです。

多数の数値の和を計算する場合、絶対値の小さいものから順に足し合わせることで、情報落ちの影響を軽減できることがあります。より洗練された手法として、後述する Kahan summation（補正加算） があります。

特殊値：NaNとInfinity

浮動小数点数には、通常の数値以外に特殊な値が存在します。これらはシミュレーション中に意図せず発生することがあり、その挙動を理解しておくことが重要です。

Infinity（無限大）

f64::INFINITY（正の無限大）と f64::NEG_INFINITY（負の無限大）は、表現可能な範囲を超えた値を示します。

• オーバーフロー（非常に大きな数の演算）で発生
• 正の数をゼロで除算した場合に発生

NaN（Not a Number）

f64::NAN は、数学的に定義できない演算結果を示します。

• 0.0 / 0.0（ゼロ除算の不定形）
• (-1.0_f64).sqrt()（負の数の平方根）
• f64::INFINITY - f64::INFINITY（無限大同士の減算）

fn main() {
    // Infinityの発生
    println!("1.0 / 0.0 = {}", 1.0_f64 / 0.0);           // inf
    println!("1.0e308 * 10.0 = {}", 1.0e308_f64 * 10.0); // inf (オーバーフロー)

    // NaNの発生
    println!("0.0 / 0.0 = {}", 0.0_f64 / 0.0);           // NaN
    println!("(-1.0).sqrt() = {}", (-1.0_f64).sqrt());   // NaN

    // NaNの特殊な性質
    let nan = f64::NAN;
    println!("\nNaN == NaN は {}", nan == nan);  // false（重要！）
    println!("NaN.is_nan() = {}", nan.is_nan()); // true

    // 有限性のチェック
    let inf = f64::INFINITY;
    println!("\n1.0.is_finite() = {}", 1.0_f64.is_finite());   // true
    println!("inf.is_finite() = {}", inf.is_finite());         // false
    println!("nan.is_finite() = {}", nan.is_finite());         // false
}

Warning

NaNの伝播に注意

NaNを含む演算の結果はすべてNaNになります。シミュレーション中に一度でもNaNが発生すると、以降のすべての計算結果が「汚染」されてしまいます。

let nan = f64::NAN;
println!("{}", nan + 1.0);   // NaN
println!("{}", nan * 0.0);   // NaN
println!("{}", nan.sin());   // NaN

結果が期待と異なる場合や、シミュレーションが発散した場合は、is_nan() や is_finite() を用いて中間値をチェックすることが有効なデバッグ手法です。

let nan = f64::NAN;
println!("{}", nan.max(0.0)); // 0（NaNではない！）
println!("{}", nan.min(0.0)); // 0（NaNではない！）

Note

オーバーフローとアンダーフロー

f64 が表現できる最大値（約 ）を超えるとInfinityになります。一方、最小の正規化数（約 ）より小さい正の値は、ゼロまたは精度が低下した「非正規化数」になります（アンダーフロー）。

非常に大きな値や小さな値を扱う計算では、対数スケールでの計算や、適切なスケーリング（無次元化）を検討してください。

誤差への対策

これらの誤差が避けられないものである以上、それらを念頭に置いたプログラミングが求められます。

浮動小数点数の比較

前述の通り、== を用いた直接比較は危険です。代わりに、2つの数値の差がごく小さい値（許容誤差、イプシロンとも呼ばれる）未満であるかどうかで判定するのが一般的です。

絶対誤差による比較

最も単純な方法は、2つの数値の差の絶対値が許容誤差未満かを確認する方法です。

fn main() {
    let a: f64 = 0.1 + 0.2;
    let b: f64 = 0.3;

    let tolerance = 1e-10; // 許容誤差

    // 絶対誤差による比較
    if (a - b).abs() < tolerance {
        println!("絶対誤差の範囲で a と b は等しい");
    } else {
        println!("絶対誤差の範囲で a と b は等しくない");
    }
}

相対誤差による比較

絶対誤差による比較は、比較する数値が非常に大きい場合や非常に小さい場合にうまく機能しないことがあります。例えば、1.0e10 と 1.0e10 + 1.0 を比較する場合、その差 1.0 は 1e-10 よりもはるかに大きいですが、値のスケールから考えれば両者は「ほぼ等しい」と見なしたいかもしれません。

このような場合には、相対誤差を用いるのが有効です。

fn relative_eq(a: f64, b: f64, tolerance: f64) -> bool {
    // 注意: NaNが渡された場合、この関数はfalseを返す
    // （NaNとの比較はすべてfalseになるため）

    // ゼロ除算を避ける
    if a == 0.0 || b == 0.0 {
        return (a - b).abs() < tolerance;
    }
    (a - b).abs() / a.abs().max(b.abs()) < tolerance
}

この方法では、差を値の大きさで正規化するため、スケールに依存しない比較が可能になります。

Note

マシンイプシロン f64::EPSILON

Rustの標準ライブラリには、f64::EPSILON という定数が用意されています。これは「1.0と、1.0より大きい次の表現可能なf64の値との差」と定義されており、マシンイプシロンと呼ばれます。 ただし、この値は非常に小さく（約 ）、比較対象の数値のスケールによっては実用的でない場合があります。通常は、問題の性質に応じて適切な許容誤差 tolerance を独自に設定する必要があります。

実践的な比較: approx クレート

多くの場合、上記のような比較関数を自前で実装するよりも、approx のような実績のあるクレートを利用するのが最善です。このクレートは、絶対誤差、相対誤差、およびULP（Unit in the Last Place）に基づいた、頑健で使いやすい比較関数を提供します。

Cargo.toml に以下を追加して approx を導入できます。

[dependencies]
approx = "0.5"

そして、assert_relative_eq! マクロなどを使ってテストを書くことができます。

use approx::assert_relative_eq;

fn main() {
    let a: f64 = 0.1 + 0.2;
    let b: f64 = 0.3;

    // a と b が相対誤差の範囲で等しいことを表明する
    assert_relative_eq!(a, b, max_relative = 1e-10);

    println!("approxクレートを使って a と b がほぼ等しいことを確認しました。");
}

シミュレーションコード内の条件分岐などで比較を行いたい場合は、relative_eq!マクロやabs_diff_eq!マクロがブール値を返すため便利です。

本書では、テストコードなどで浮動小数点数の比較を行う際に approx クレートを積極的に利用します。

補正加算（Kahan Summation）

多数の数値を足し合わせる際、単純に加算を繰り返すと情報落ちによる誤差が蓄積します。Kahan summation（補正加算、カハンの加算アルゴリズム）は、誤差を追跡・補正しながら加算を行うことで、この問題を軽減するアルゴリズムです。

/// Kahan summationによる総和計算
fn kahan_sum(values: &[f64]) -> f64 {
    let mut sum = 0.0;
    let mut compensation = 0.0; // 累積誤差の補正項

    for &value in values {
        let y = value - compensation;  // 補正を適用した値
        let t = sum + y;               // 一時的な合計
        compensation = (t - sum) - y;  // 丸め誤差を記録
        sum = t;
    }
    sum
}

fn main() {
    // 小さな値を大量に足し合わせる例
    let n = 1_000_000;
    let small_value = 0.1;
    let values: Vec<f64> = vec![small_value; n];

    // 単純な加算
    let naive_sum: f64 = values.iter().sum();

    // Kahan summation
    let kahan = kahan_sum(&values);

    // 理論値
    let expected = small_value * n as f64;

    println!("理論値:           {:.17}", expected);
    println!("単純な加算:       {:.17}", naive_sum);
    println!("Kahan summation: {:.17}", kahan);
    println!("\n単純な加算の誤差: {:.2e}", (naive_sum - expected).abs());
    println!("Kahanの誤差:      {:.2e}", (kahan - expected).abs());
}

実行結果（環境により若干異なる場合があります）：

理論値:           100000.00000000000000000
単純な加算:       100000.00000133288267534
Kahan summation: 100000.00000000000000000

単純な加算の誤差: 1.33e-6
Kahanの誤差:      0.00e0

Kahan summation は、数値積分やモンテカルロ法など、多数の値を累積する計算で特に有効です。本書の後の章でも、精度が重要な場面でこのアルゴリズムを活用します。

数値的に安定なアルゴリズムの選択

同じ問題を解くアルゴリズムでも、計算手順によって誤差の蓄積度合いが異なる場合があります。桁落ちの例で見たように、数式を適切に変形したり、より数値的に安定したアルゴリズムを選択したりすることが重要です。

高精度計算の利用

f64 の精度でも不十分な場合は、より多くのビットを使って数値を表現する「多倍長浮動小数点数」ライブラリを利用するという選択肢もあります。

これについては、付録E「高精度計算」で詳しく解説します。

まとめ

本節では、計算物理学における基本でありながら、非常に重要な浮動小数点数の誤差について学びました。

• 丸め誤差・桁落ち・情報落ちといった誤差は、浮動小数点数演算において本質的に避けられない。
• NaNやInfinityといった特殊値の発生に注意し、必要に応じて is_nan() や is_finite() でチェックする。
• == による浮動小数点数の直接比較は危険であり、許容誤差を用いた比較が必要である。
• Kahan summation などの手法を用いることで、多数の値の加算における誤差を軽減できる。
• 計算順序の工夫や、数値的に安定なアルゴリズムの選択によって、誤差の影響を軽減できる。

シミュレーションを行う際は、常にこれらの誤差が結果に与える影響を意識することが、科学的に信頼できる結論を導くための第一歩となります。