Juliaで丸めモードを指定して浮動小数点数の計算をする(したい)

2020/4/14追記 丸めモードについては、デフォルトの丸めモードで上付き丸め、下付き丸めをエミュレートする(Julia) も参照してください。 IntervalArithmetic.jl のデフォルト丸めモードは、#370 で FastRounding.jl を利用しなくなった (v0.17.0以降) ので、 IntervalArithmetic.jl に関して以下に書いてある情報は古いです。

最近、 数値計算を 区間演算 を使って誤差の評価を伴って実施する 精度保証付き数値計算 に興味が出てきて、 Warwich Tucker の Validated Numeric という本を購入し、少しつづ読んでいます。

浮動小数点数のシステムなど初歩的なことから説明されていて私にとっては非常にありがたく、異常な計算例も多く載せられていて楽しく読めています。 実際に手を動かして理解するために、本文内の MATLAB で書かれた区間演算のプログラムを Julia で実装しようとしていますが、 丸めモードを指定した計算を行うところでつまづいてしまいました。

今回は、頭の整理を兼ね、具体的な例からスタートし、なぜ丸めモードを指定した計算が必要となるか説明するとともに、 Julia における丸めの制御について調べた内容を書いておきたいと思います。 私はこの分野の専門家ではないので、間違っている場所があれば教えていただけると助かります。

使用している Julia のバージョンは 1.3.0 です。将来的に丸めの扱いが変更される可能性があるので注意してください。

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.3.0 (2019-11-26)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

0.1 * 3.0 を計算してみる

$\mathbb{F}^*$ を浮動小数点全体の集合 (今回はFloat64 としておきます)、 $\mathbb{F}$ を$\mathbb{F}^*$から -Inf, +Inf を除いた有限な浮動小数点全体の集合、 $\mathbb{R}$ を実数全体の集合、$\mathbb{R}^*$ を $\mathbb{R}$ に正と負の無限大を付け加えた集合とします (無限大が関係してくる細かい説明は今回は省略します)。 $\mathbb{F}$ は $\mathbb{R}$、 $\mathbb{F}^*$ は $\mathbb{R}^*$ の有限部分集合です。

例として、 $0.1 \in \mathbb{R}$ と $3.0 \in \mathbb{R}$ の積を (浮動小数点上で) 評価したいとします。 言うまでもなく真の値は $0.3 \in \mathbb{R}$ です。 一方、REPL 上で 0.1 * 3.0 と 0.3 を比較すると、

julia> 0.1 * 3.0 == 0.3
false

となりイコールにはなりません。@printf マクロを使って小数点以下17桁まで表示すると、確かに違います。

julia> using Printf

julia> @printf("%17.17f", 0.1 * 3.0)
0.30000000000000004
julia> @printf("%17.17f", 0.3)
0.29999999999999999

Base.bitstring を使ってビット表現を確かめても、

julia> bitstring(0.1 * 3.0)
"0011111111010011001100110011001100110011001100110011001100110100"

julia> bitstring(0.3)
"0011111111010011001100110011001100110011001100110011001100110011"

最後の3桁が違いますね。なぜこのようなことが起こるのでしょうか？

丸め

証明は省略しますが、$1/3$ が有限桁の10進数小数で表現できないのと同様、 $0.1$ は有限桁の(2進)浮動小数点数では表せません。$0.3$ も同様です。

そのため、$0.1 \in \mathbb{R}$ はコンピューターで扱うために $\mathbb{F}$ の元で近似されることになります。 この近似の方法を定めるのが丸め (Rounding) で、書き方は違いますが本の中で次のように定義されています。

定義 「丸め」 とは、写像 $\bigcirc : \mathbb{R}^* \rightarrow \mathbb{F}^*$ で、
(R1) $\bigcirc$ は $\mathbb{F}^*$ 上では恒等写像。つまり、任意の $x \in \mathbb{F}^*$ に対して $\bigcirc (x) = x$
(R2) $\bigcirc$ は大小関係を保存する。つまり、$x, y \in \mathbb{R}^*$ が $x \le y$ を満たすなら $\bigcirc(x) \le \bigcirc(y)$ が成立する
の二つの条件を満たすものである。

一つ目の条件は $\mathbb{F}^*$ でもともと表せる数に関しては余計な操作を行わないことを意味し、 二つ目の条件は近似した時に大小関係が入れかわらないことを保証するもので、どちらも非常に自然な要請です。

丸めにはいくつか種類がありここでは全て紹介できませんが、いくつか基本的な丸めを挙げます。

はそれぞれ下付き丸め (Rounded down), 上付き丸め (Rounded up) と呼ばれます。

これらは勿論先ほどの丸めの定義を満たし、

が常に成立しています。これと条件 (R1), (R2) を組み合わせることにより、任意の丸め $\bigcirc$ に対して

が成立することがわかります。また、$x \in \mathbb{F}^*$ であれば、 $\bigcirc (x) = x = \triangledown (x) = \triangle(x)$ ですし、 $x \notin \mathbb{F}^*$ であれば、$\triangledown(x) < x < \triangle (x)$ が成り立つとともに $\bigcirc (x) = \triangledown(x)$ または $\bigcirc (x) = \triangle(x)$ となることも示せます。 つまり、丸めとは、浮動小数点数で表せないそれぞれの数に対して切り上げるか切り下げるかを条件 (R2) を守るように定めることと同じです。

浮動小数点数で表現できない数に対する Julia のデフォルト丸めモードは RoundNearest で、 最も近い浮動小数に丸めると言うものです。 ちょうど中間になった時が気になりますが、切り上げた時と切り下げた時のうち、最後の桁が偶数 (つまり0) になるものを選びます。 そのためこの丸めモードはNearest Evenと呼ばれることもあります。 なぜ偶数に丸めるかについての説明は省略しますが本には載っているので気になる方は参照してください。

今後、このデフォルトの丸めを $\square$ で表すことにします。0.1, 3.0, 0.3 は、 実態はそれぞれ $\square(0.1), \square(3.0) = 3.0, \square(0.3) \in \mathbb{F}^*$ だったということになります。

浮動小数点数同士の演算

次に $\mathbb{R}^*$ 上の二項演算 $+, -, \times, /, \sqrt{}$ とそれに対応する $\mathbb{F}^*$ 上の二項演算 $\oplus, \ominus, \otimes, \oslash, \odot$ の違いを考えます。 本当は四角の囲み文字を使いたかったのですが、ご勘弁ください。ルートの丸文字も出せませんでした。

問題は $+, -, \times, /, \sqrt{}$ を $\mathbb{F}^*$ に制限しても結果が $\mathbb{F}^*$ になるとは限らないと言うことです。 例えば、$2^{100}, 1 \in \mathbb{F}^*$ ですが、$2^{100} + 1 \notin \mathbb{F}^*$ です。 Julia には直前、直後の浮動小数点数を計算する関数 prevfloat, nextfloat が存在しますので $2^{100}$の次の数を計算してみましょう。

julia> x = 2.0^100
1.2676506002282294e30

julia> @printf("%17.17f", x)
1267650600228229401496703205376.00000000000000000
julia> @printf("%17.17f", nextfloat(x))
1267650600228229682971679916032.00000000000000000

なので、$2^{100} + 1$ は確かにこの二つの隙間に落ちてしまっているので表現不能です。計算すると、

julia> x + 1.0
1.2676506002282294e30

julia> x + 1.0 == x
true

となって元々の値と同じになっています。

つまりは計算した結果も丸められているわけですが、ここで丸める前と丸めた後の関係が気になります。 Juliaの Rounding modes を見ると丸めモードは IEEE754 standard に従っていると言うことなので、 規格を本来は参照すべきだと思いますが、ノンメンバーだとPDFを購入するのに$100かかるようなので、断念しました。

書かれたのは1997年と古いようですが、IEEE 754の策定に携わったKahanのレクチャーノート Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic を読むと、四則演算や平方根丸め無しの場合に計算した結果を Nearest Even で丸めるという内容のことが書かれています。 別の本「精度保証付き数値計算の基礎」（これも勢いで購入）にも同様のことが書いてあった (p.p. 17) ので、信じることにします。

よって、

の時に

ということになります。

長々と書いてしまいましたが、0.1 * 3.0 $= \square(\square(0.1) \times \square(3.0))$, 0.3 $= \square(0.3)$ なので、0.1 * 3.0 と 0.3 は近い値にはなるものの、$\square(0.1) \times \square(3.0) \neq 0.3$ でそれぞれ別の数に丸められてしまうということでした。

区間演算

次に、同様の計算を区間計算で (正確な定義はスキップして) 行ってみます。 まず、$0.1$ を含む、端点が共に $\mathbb{F}$ の元で幅がなるべく小さい区間を探します。

0.1 と書いた時点で浮動小数点数に変換されてしまうので、別の方法で $0.1$ を表現する必要があります。 入力を文字列 "0.1" にしてパースする方法もありますが、 幸いなことに有理数型と丸め方向つきの浮動小数点数コンストラクタが利用できます。

julia> x_down, x_up = Float64(1//10, RoundDown), Float64(1//10, RoundUp);

julia> @printf("[%17.17f, %17.17f]", x_down, x_up)
[0.09999999999999999, 0.10000000000000001]

$3.0$ は浮動小数点数として誤差なく表現できるので、幅0 (thin) の区間として表わせます。

julia> Float64(3//1, RoundDown) == Float64(3//1, RoundUp)
true

こうして $0.1$ と $3.0$ がそれぞれ含まれる区間 [x_down, x_up], [3.0, 3.0] が得られたので答えが取り得る範囲を考えて、真の数が含まれる区間の計算結果を [x_down * 3.0, x_up * 3.0] ……とやると(一般的には)まずいです。

今回はたまたま

julia> @printf("[%17.17f, %17.17f]", x_down * 3.0, x_up * 3.0)
[0.29999999999999999, 0.30000000000000004]

となってうまくいきましたが、例えば、 $1.0 \in \mathbb{F}, 2^{-54} \in \mathbb{F}$ に対し、 $1.0 + 2^{-54}$ を意図して同様の区間演算 [1.0, 1.0] + [2^-54, 2^-54] を先ほどと同様に単純に行うと、

julia> @printf("%17.17f", 1.0 + 2^-54)
1.00000000000000000
julia> 1.0 + 2^-54 == 1.0
true

と $2^{-54}$ が $1$ に対して小さすぎて桁落ちが発生し、 得られる区間は [1.0, 1.0] となってしまい正しい答えである $1.0 + 2^{-54}$ は含まれません。

つまり、今までに定義した記号を使って表すと、

と定義してしまってはダメで、

と計算して答えが確実に含まれるようにする必要があります。ここで必要となるのは、デフォルトの丸めモードとは異なる丸めモードを用いた加算です。

非常に長くなってしまいましたが、この計算をやりたいがために丸めモードを変えたいわけです。

Julia における方向丸め付き演算のサポート

「さあ、いよいよ Julia で方向丸め付き演算をやっていくか、C++だと丸めモードの変更は面倒な感じがする が、 Julia には Base.Rounding.setrounding という関数が用意されていて楽すぎる。Julia最強！」

……とはならないのが人生の辛い所です。

setrounding の説明をよく読むと、T == BigFloat しかサポートされていないと書かれています。 実は昔は Float32, Float64 に対しても使えていたようなのですが、 色々と問題があり(最後のリンク参照)、 experimentalと言う但し書きがドキュメントに書かれるようになり、 最終的に廃止されたようです。

そのためブログの記事や Twitter の投稿で「Juliaでは浮動小数点の方向丸め付き演算がサポートされている」と書かれているものは、おそらく、 Float32, Float64 に対する setrounding が使えるようになっていた時のものでしょう。

サポートされていない理由については、どうもLLVMで浮動小数点の丸めの制御ができないため、と言うことらしいです(未確認)。

Julia で利用できるパッケージ

一方、Juliaには区間演算のライブラリ IntervalArithmetic.jl があり、頻繁にアップデートされています。 丸めは FastRounding.jl の機能を使っているようでした。

他にも DirectedRoundings.jl や SetRounding.jl を発見しましたが、あまりメンテされていないようでした。

最近点丸めによる方向付き丸め演算のエミュレート

言語でサポートされていない方向付き丸めがパッケージを使えばできるということで意味がわからなくなってしまったのですが、 FastRounding.jl は、ErrorfreeArithmetic.jl で実装されているエラーフリー変換と言うものを用いて、 デフォルトの丸めモードの演算のみを利用して方向付き丸めをエミュレートしているようです。

エミュレートの基本原理に関しては、 最近点丸めによる方向付き丸めのエミュレート に詳細に書いてあるのですが、 ここでは簡単に加算に関して FastRounding.jl が行っているエラーフリー変換とエミュレートについて説明してみたいと思います。

加算のエラーフリー変換は、別名 twosum とも呼ばれるアルゴリズムで、$a, b \in \mathbb{F}$ に対して, ペア $(x, y) = \text{twosum}(a, b)$ を返し、

が厳密に成立するとされるものです。

Juliaで書くとほぼ「最近点丸めによる方向付き丸めのエミュレート」に書かれている Python のコードと変わらず、

function twosum(a, b)
    x = a + b
    tmp = x - a
    y = (a - (x - tmp)) + (b - tmp)
    x, y
end

となります。 私が確認したタイミングでは、ErrorfreeArithmetic.jl の two_sum の定義も同様になっていました。

今まで出てきた例の場合も含めて何通りか例を計算してみると

julia> twosum(0.1, 0.2)
(0.30000000000000004, -2.7755575615628914e-17)

julia> 2.0^100
1.2676506002282294e30

julia> twosum(2.0^100, 1.0)
(1.2676506002282294e30, 1.0)

julia> 2^-54
5.551115123125783e-17

julia> twosum(1.0, 2^-54)
(1.0, 5.551115123125783e-17)

となって正しく分解されています。

$a, b$ の絶対値について条件を課したアルゴリズム fasttwosum について、 オーバーフローが発生しない場合に等式が成立することが「精度保証付き数値計算の基礎」に定理2.5として証明されていました。

twosum を用いた方向付き丸め演算のエミューレートは、上付き丸めと下付き丸めを

function add_up(a, b)
    x, y = twosum(a, b)
    ifelse(y > 0, nextfloat(x), x)
end

function add_down(a, b)
    x, y = twosum(a, b)
    ifelse(y < 0, prevfloat(x), x)
end

と定義するもので、FastRounding.jlの対応箇所はここです。

julia> add_up(2.0^100, 1.0)
1.2676506002282297e30

julia> add_down(2.0^100, 1.0)
1.2676506002282294e30

julia> add_up(1.0, 2^-54)
1.0000000000000002

julia> add_down(1.0, 2^-54)
1.0

見た感じ、良さそうです。FastRounding.jl の関数も呼び出してみます。

julia> using FastRounding

julia> add_round(2.0^100, 1.0, RoundUp)
1.2676506002282297e30

julia> add_round(2.0^100, 1.0, RoundDown)
1.2676506002282294e30

julia> add_round(1.0, 2^-54, RoundUp)
1.0000000000000002

julia> add_round(1.0, 2^-54, RoundDown)
1.0

しかし、これで全て安心とはならず、救われない命があるようなのです……

エミュレーションにおける問題

Error Free Transformation (EFT) is NOT error-free という記事に、 エラーフリー変換においてオーバーフローやアンダーフローが発生する例や、twosum の乗算バージョンの twoproduct に関して誤差が生じる例が書かれていました。

気になったので、先ほどの twosum に対して記事に書かれている例を試してみると

julia> a = 3.5630624444874539e+307
3.563062444487454e307

julia> b = -1.7976931348623157e+308
-1.7976931348623157e308

julia> twosum(a, b)
(-1.4413868904135704e308, NaN)

となってやはり正しく計算されませんでした。

「最近点丸めによる方向付き丸めのエミュレート」にはこれを回避する方法が書いてあったので、実装してみます。

function new_twosum(a, b)
    x = a + b
    if abs(a) > abs(b)
        tmp = x - a
        return x, b - tmp
    else
        tmp = x - b
        return x, a - tmp
    end
end

function new_add_up(a, b)
    x, y = new_twosum(a, b)
    ifelse(y > 0, nextfloat(x), x)
end

function new_add_down(a, b)
    x, y = new_twosum(a, b)
    ifelse(y < 0, prevfloat(x), x)
end

実行結果は、もちろん今まで計算できていたものは同じ結果になり、

julia> new_twosum(0.1, 0.2)
(0.30000000000000004, -2.7755575615628914e-17)

julia> new_twosum(2.0^100, 1.0)
(1.2676506002282294e30, 1.0)

julia> new_twosum(1.0, 2^-54)
(1.0, 5.551115123125783e-17)

julia> new_add_up(2.0^100, 1.0)
1.2676506002282297e30

julia> new_add_down(2.0^100, 1.0)
1.2676506002282294e30

julia> new_add_up(1.0, 2^-54)
1.0000000000000002

julia> new_add_down(1.0, 2^-54)
1.0

先ほど計算できなかった例が、

julia> new_twosum(a, b)
(-1.4413868904135704e308, 9.9792015476736e291)

julia> new_add_up(a, b)
-1.4413868904135702e308

julia> new_add_down(a, b)
-1.4413868904135704e308

と計算できるようになりました。

しかしながら ErrorfreeArithmetic.jl の実装も残念ながらオーバーフローが考慮されていないので、

julia> using ErrorfreeArithmetic

julia> two_sum(a, b)
(-1.4413868904135704e308, NaN)

となり、従ってそれを利用した丸めも

julia> using FastRounding

julia> add_round(a, b, RoundUp)
-1.4413868904135704e308

julia> add_round(a, b, RoundDown)
-1.4413868904135704e308

と丸めが正しく行われていません。記事に載っていた例については加算しか確認していませんが、他にも似たような現象が発生している可能性があるので利用には注意が必要だと思います。 (ライブラリの作者の方々には申し訳ないですが)

基本的な演算に関してのエミュレートを行うだけでもかなり気を使う必要があり、漏れなく実装するのは容易ではないようです。

しかしながら、デフォルトの丸めモードだけを使って方向付き丸めのエミュレートを行う発想自体は非常に面白く感じたので、時間のある時に他の演算に関しても実装してみたいです。 (パッケージとして通用するレベルで厳密な実装やテストを書くのは大変そうですが……)

結論(2020/3時点)

Julia で方向付き丸めを指定した浮動小数点数計算はコア言語としては廃止されているが、エミュレートにより丸めを行うパッケージ FastRounding.jl が存在します。

しかしながら FastRounding.jl は極端な入力を与えた場合に正しい答えを返さない可能性があり、利用の際には注意が必要、と言う結論になりました。

長くなりましたが最後に参考資料を載せておきます。

参考資料、リンク

精度保証付数値計算

• Tucker, W. Validated Numerics: A Short Introduction to Rigorous Computations, Princeton University Press, 2011
• 大石 新一. 精度保証付き数値計算の基礎, コロナ社, 2018
• Wikipedia. 精度保証付き数値計算
• JuliaIntervals / ValidatedNumerics.jl
• JuliaIntervals / ValidatedNumerics.jl / Julia rounding problem #165

IEEE 754

• Kahan, W. Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic
• Wikipedia. IEEE 754
• ”Validated Numerics" の 1.5章
• IEEE Computer Society. 754-2019 - IEEE Standard for Floating-Point Arithmetic

Julia言語の丸めモード

• Integers and Floating-Point Numbers > Rounding modes
• Base.Rounding.setrounding
• Base.Rounding.RoundingMode
• JuliaLang / julia / setrounding regression for 0.5.0-rc1 #17926
• JuliaLang / julia / what to do with setrounding? #26935
• JuliaLang / julia / Deprecate setrounding #27166
• o108minmin. JuliaのFloatとRoundingの奇妙な挙動[追記]

方向付き丸めのエミュレート, Error Free Transformation

• 柏木 雅英. 最近点丸めによる方向付き丸めのエミュレート
• 柏木 雅英. 最近点丸めのみによる方向付き丸めのエミュレート
• Takeshi O, Siegfried R, and Shin’ichi O, Accurate Sum and Dot Product, SIAM J. Sci. Comput., 26(6), 1955–1988 (http://www.oishi.info.waseda.ac.jp/%7Eoishi/papers/OgRuOi05.pdf でプレプリントが読める)
• 「精度保証付き数値計算」の 2.2 章
• 柏木 雅英. Error Free Transformation (EFT) is NOT error-free
• JeffreySarnoff / FastRounding.jl
• JeffreySarnoff / ErrorfreeArithmetic.jl

方向つき丸めのJuliaパッケージ(未調査)

• JeffreySarnoff / DirectedRoundings.jl
• JuliaIntervals / SetRounding.jl

区間演算

• Validated Numerics の 2章
• 柏木 雅英. 区間演算の実装について(1)
• 柏木 雅英. 区間演算の実装について(2)
• @mod_poppo. 週刊 代数的実数を作る #5 区間演算と計算可能実数
• JuliaIntervals / IntervalArithmetic.jl
• JuliaIntervals / IntervalArithmetic.jl / Bugs with tiny numbers #215