Rust で quad-double の加算を試す

Tech
Share

最近 Rust という言語を学んでいます。システムプログラム向けということらしいですが、私はパフォーマンスからして数値計算にも向いているのではないかと勝手に考えています。

行列計算プログラムでも書けばいいかもしれませんが、偶然出会った double-double 演算に興味を持ったので、 Rust で quad-double の加算を実装してきました。 C++ による実装はすでにありますので参考にしつつ。

参考論文はこちら、実装はこちら。

quad-double 演算について

quad-double （以下 QD と書く）演算は、倍精度浮動小数点（ double ）を 4 つ使うことにより、精度の高い演算を実現します。 2 つ使ったものが double-double 演算です。 4倍精度浮動小数点（ gcc での long double ）や 8 倍精度浮動小数点数とは違って値のとれる範囲は増えませんが、精度の必要な計算に対するソフトウェア的な手法として使われることがあるようです。 Rust では double は f64 に当たりますが、たまに Rust への要望に出てくる f128 は 4 倍精度のほうだと思われます。

データ構造の実装

QD 演算の実装に際し、論文に対応するよう 3 つの構造体と関数を作成しました。

1
2
3
4
5
6
7
8
9
10
11
#[derive(Debug, Copy, Clone)]
struct SumError {
    sum: f64,
    error: f64,
}

#[derive(Debug, Copy, Clone)]
struct ExtendedQuadDouble(f64, f64, f64, f64, f64);

#[derive(Debug, Copy, Clone)]
struct NormalizedQuadDouble(f64, f64, f64, f64);
0
// omit
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#[inline]
fn renormalize(a: ExtendedQuadDouble) -> NormalizedQuadDouble {
    let SumError { sum: s, error: t4 } = quick_two_sum(a.3, a.4);
    let SumError { sum: s, error: t3 } = quick_two_sum(a.2, s);
    let SumError { sum: s, error: t2 } = quick_two_sum(a.1, s);
    let SumError { sum: t0, error: t1 } = quick_two_sum(a.0, s);
    let t = [t0, t1, t2, t3, t4];

    let mut b = [t[0], t[1], t[2], t[3]];
    let mut s = t[0];
    let mut k = 0;
    for i in 1..5 {
        let r = quick_two_sum(s, t[i]);
        s = r.sum;
        let e = r.error;
        if e != 0.0 {
            b[k] = s;
            s = e;
            k = k + 1;
        }
    }
    NormalizedQuadDouble(b[0], b[1], b[2], b[3])
}

SumError は f64 の加算の結果と情報落ちのペアです。 QD 数は演算直後は 5 つの f64 で表される展開形式 ExtendedQuadDouble で表され、それを唯一に定まる正規化表現 NormalizedQuadDouble に変換します。

double 同士の加算

いつもの double 同士の加算では暗黙の了解で誤差を無視してしまっていますが、誤差について考えると、条件によって誤差の算出に必要な演算回数が違う 2 つの加算が生まれます。

13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#[inline]
fn quick_two_sum(a: f64, b: f64) -> SumError {
    let s = a + b;
    let e = b - (s - a);

    SumError { sum: s, error: e }
}

#[inline]
fn two_sum(a: f64, b: f64) -> SumError {
    let s = a + b;
    let v = s - a;
    let e = (a - (s - v)) + (b - v);

    SumError { sum: s, error: e }
}

quick_two_sum() は ∣a∣≥∣b∣|a|\geq|b| が成り立つとき、 two_sum() はいつでも使えます。正規化関数 renormalize() では quick_two_sum()を使うことで演算回数を減らしています。

QD 数との加算

ここでは QD 数と double の加算と QD 数同士の加算を実装しています。

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#[inline]
fn qd_d_add(a: NormalizedQuadDouble, b: f64) -> NormalizedQuadDouble {
    let SumError { sum: t0, error: e } = two_sum(a.0, b);
    let SumError { sum: t1, error: e } = two_sum(a.1, e);
    let SumError { sum: t2, error: e } = two_sum(a.2, e);
    let SumError { sum: t3, error: t4 } = two_sum(a.3, e);

    renormalize(ExtendedQuadDouble(t0, t1, t2, t3, t4))
}

#[inline]
fn double_accumulate(u: f64, v: f64, x: f64) -> (f64, f64, f64) {
    let SumError { sum: s, error: mut v } = two_sum(v, x);
    let SumError { sum: mut s, error: mut u } = two_sum(u, s);
    if u == 0.0 {
        u = s;
        s = 0.0;
    }
    if v == 0.0 {
        v = u;
        u = s;
        s = 0.0;
    }

    (s, u, v)
}

#[inline]
fn qd_add_accurate(a: NormalizedQuadDouble, b: NormalizedQuadDouble) -> NormalizedQuadDouble {
    // merge sort
    let a = [a.0, a.1, a.2, a.3];
    let b = [b.0, b.1, b.2, b.3];
    let mut x = [0.0; 8];

    let mut i = 0; // for a
    let mut j = 0; // for b
    for ptr in 0..8 {
        if b.len() <= j || a.len() > i && a[i].abs() > b[j].abs() {
            x[ptr] = a[i];
            i += 1;
        } else {
            x[ptr] = b[j];
            j += 1;
        }
    }
    let x = x;

    let mut u = 0.0;
    let mut v = 0.0;
    let mut k = 0;
    let mut i = 0;
    let mut c = [0.0; 4];
    while k < 4 && i < 8 {
        let r = double_accumulate(u, v, x[i]);
        let s = r.0;
        u = r.1;
        v = r.2;

        if s != 0.0 {
            c[k] = s;
            k = k + 1;
        }
        i = i + 1;
    }
    if k < 2 {
        c[k + 1] = v;
    }
    if k < 3 {
        c[k] = u;
    }

    renormalize(ExtendedQuadDouble(c[0], c[1], c[2], c[3], 0.0))
}

実行例

128
129
130
131
132
133
134
135
136
137
138
139
140
fn main() {
    // lose data
    println!("{:.10}", 3E-9_f64 + 7E+9_f64);

    let qd = NormalizedQuadDouble(3E-9_f64, 0.0, 0.0, 0.0);
    let qd_result_1 = qd_d_add(qd, 7E+9_f64);
    println!("{:?}", qd_result_1);

    let qd_big = NormalizedQuadDouble(3E+18_f64, 5E+9_f64, 0.0, 0.0);
    let qd_small = NormalizedQuadDouble(7E-9_f64, 9E-18_f64, 0.0, 0.0);
    let qd_result_2 = qd_add_accurate(qd_big, qd_small);
    println!("{:?}", qd_result_2);
}
1
2
3
4
7000000000.0000000000
NormalizedQuadDouble(7000000000, 0.000000003, 0, 0)
NormalizedQuadDouble(3000000005000000000, 0.0000000070000000089999995, 0.00000000000000000000000031145
969116572085, 0)

まず情報落ちの例として、 3E-9_f64 + 7E+9_f64 の結果が 7000000000.0000000000 となり、 double 同士の加算では 3E-9_f64 の情報が消えていることがわかります。 QD 数と double の加算については、しっかり 3E-9_f64 の情報が保存されていることがわかります。 QD 数同士の加算については、多少の誤差が出てしまいましたが、なんとか 3E+18_f64 、 5E+9_f64 、 7E-9_f64 、 9E-18_f64 の存在がわかります。