matsueushi | デフォルトの丸めモードで上付き丸め、下付き丸めをエミュレートする(Julia)

Juliaで丸めモードを指定して浮動小数点数の計算をする(したい) で色々とJuliaの丸めモードについて調べていましたが、今回はその続きです。

やはりどうしてもJuliaの Float32 や Float64 に対して上付き丸め、下付き丸めを可能な限り正確に安定して行いたくなったので、「最近点丸めによる方向付き丸めのエミュレート」を参考に、デフォルトの丸めモードのみを使って上付き丸め、下付き丸めをエミュレートするJuliaのパッケージを作成しました。

RoundingEmulator.jl

https://github.com/matsueushi/RoundingEmulator.jl

Registratorにも登録したので、

] add RoundingEmulator

でインストールして使えます。

定義されているのは add, sub, mul, div, sqrt に対して上付き丸め _up と下付き丸め _down と2つと非常なシンプルなものです。

julia> using RoundingEmulator

julia> add_up(0.1, 0.2)
0.30000000000000004

julia> add_down(0.1, 0.2)
0.3

julia> sub_up(-0.1, 0.2)
-0.3

julia> sub_down(-0.1, 0.2)
-0.30000000000000004

julia> mul_up(0.1, 0.2)
0.020000000000000004

julia> mul_down(0.1, 0.2)
0.02

julia> div_up(1.0, 3.0)
0.33333333333333337

julia> div_up(1.0, 3.0)
0.33333333333333337

julia> sqrt_up(2.0)
1.4142135623730951

julia> sqrt_down(2.0)
1.414213562373095

https://github.com/JeffreySarnoff/FastRounding.jl/blob/03ff386d36aa7de101f22ca740748a31e57942c0/src/FastRounding.jl#L187-L194 の例も正しく計算できていました。

julia> sqrt_up(3.9036066558023176e-154)
1.9757547053726885e-77

julia> sqrt_down(3.9036066558023176e-154)
1.975754705372688e-77

詳細はレポジトリやドキュメンテーションを参照して下さい。 Issueやプルリクは大歓迎です。

基本的には「最近点丸めのみによる方向付き丸めのエミュレート」に沿った実装なので、エミュレートのロジックについてはそちらを参照してください。

これで、上付き丸めや下付き丸めができるようになったので、当初やろうと思っていた区間演算にもチャレンジしてみたいです。

以下は作成していて気づいた点です。

`add12`, `mul12`

エクスポートされていない関数ですが、Julia の base/twiceprecision.jl には加算、乗算のエラーフリー変換の関数 add12 と mul12 が定義されています。前回言及したオーバーフロー対策もなされていました。

julia> using Base: add12, mul12

julia> a = 3.5630624444874539e+307
3.563062444487454e307

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

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

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

julia> a = 6.929001713869936e+236
6.929001713869936e236

julia> b = 2.5944475251952003e+71
2.5944475251952003e71

julia> mul12(a, b)
(1.7976931348623157e308, -1.0027614963959625e291)

julia> mul12(b, a)
(1.7976931348623157e308, -1.0027614963959625e291)

`setrounding_raw`

deprecatedになった setrounding は、強引に

julia> using Base.Rounding: setrounding_raw, to_fenv

julia> 0.1 + 0.2
0.30000000000000004

julia> setrounding_raw(Float64, to_fenv(RoundDown))
0

julia> 0.1 + 0.2
0.3

julia> setrounding_raw(Float64, to_fenv(RoundUp))
0

julia> 0.1 + 0.2
0.30000000000000004

julia> setrounding_raw(Float64, to_fenv(RoundNearest))
0

と呼び出せば使えて、確実に動くことが保証されていないですが、テストを行うにあたっては、Windows, MacOS, LinuxのCIで自分で計算した値と、setrounding_raw で丸めモードを変更した時の値が等しくなることを特殊な値 (zero(T), floatmax(T), flaotmin(T) など)とランダムにサンプリングした値に対して確かめる、という方法を現在は取っています。

確実に正しいことがわかっている入力と出力のペアからテストを作成するのが一番良いと思うのですが、今回はそこまでやっていません。

Signed zero

符号付きゼロもの符号も含めて正確に再現する場合、下付き丸めによる加算の定義で注意する必要があります。上付き丸めによる加算は

function add_up(a, b)
    x, y = add12(a, b)
    if isinf(x)
        ifelse(x == typemin(x) && isfinite(a) && isfinite(b), -floatmax(x), x)
    else
        y > zero(y) ? nextfloat(x) : x
    end
end

と定義できますが、下付きの丸めに関しては、0どうしの計算を行った際に符号を変える必要があります。これは和の合計が0になる計算 (exact zero sum) の挙動が下付き丸めとそれ以外で異なるためです。 IEEE 754の6.3章に下付き丸め(roundTowardNegative)の時のみ、exact zero sumは-0になると書いてあります。ただ、\(x + x\) の符号は \(x\) の符号と同じになると定められているので、表にすると、

Rounding	\((+0)+(+0)\)	\((-0)+(+0)\)	\((+0)+(-0)\)	\((-0)+(-0)\)
roundTiesToEven (`RoundNearest`)	\(+0\)	\(+0\)	\(+0\)	\(-0\)
roundTowardPositive (`RoundUp`)	\(+0\)	\(+0\)	\(+0\)	\(-0\)
roundTowardNegative (`RoundDown`)	\(+0\)	\(-0\)	\(-0\)	\(-0\)

このようになります。よって、下のように修正することでゼロの符号を調整しています。

function add_down(a, b)
    x, y = add12(a, b)
    if isinf(x)
        ifelse(x == typemax(x) && isfinite(a) && isfinite(b), floatmax(x), x)
    elseif y < zero(y)
        prevfloat(x)
    else
        ifelse(iszero(x) && (signbit(a) || signbit(b)), -zero(x), x)
    end
end

計算結果:

julia> add_up(-0.0, 0.0)
0.0

julia> add_down(-0.0, 0.0)
-0.0

julia> add_up(0.0, 0.0)
0.0

julia> add_down(0.0, 0.0)
0.0

julia> add_up(-0.0, -0.0)
-0.0

julia> add_down(-0.0, -0.0)
-0.0

浮動小数点数をランダムにサンプリングする

全ての浮動小数点数から一様にサンプリングしたい時、rand を普通に呼び出すと \([0, 1)\) から一様にサンプリングされます。

https://docs.julialang.org/en/v1/stdlib/Random/index.html#Random-Numbers-1

julia> rand(Float64, 10)
10-element Array{Float64,1}:
 0.6021596191997549
 0.8268352152178551
 0.22765811638234856
 0.3813150107932535
 0.49112818878735265
 0.6424070543287013
 0.27970019663531676
 0.7316980433063847
 0.5721143543965341
 0.24735574564535145

Int64 用の rand を使って乱数を発生させ、reinterpret で Float64 に変換すれば良いです。

julia> x = rand(Int64)
7738226609100433883

julia> bitstring(x)
"0110101101100011101101000001101000100010111010000010100111011011"

julia> y = reinterpret(Float64, x)
2.0242815034259582e209

julia> bitstring(y)
"0110101101100011101101000001101000100010111010000010100111011011"

細かいことを言うと NaN, Inf なども含まれてしまいますが、今回は特別扱いはせずそのまま使っています。

`==` vs `isequal`

Julia でユニットテストを書くとき、

@test f(x) == 1.0

のように書くことが多いのですが、浮動小数点数のテストの際、NaN や符号付きゼロなどをテストしたい場合は注意が必要です。まず、

julia> NaN == NaN
false

です。また、IEEE 754で定められている符号付きゼロ -0.0 と 0.0 は

julia> bitstring(0.0)
"0000000000000000000000000000000000000000000000000000000000000000"

julia> bitstring(-0.0)
"1000000000000000000000000000000000000000000000000000000000000000"

と異なるビット列が対応していますが、== では

julia> 0.0 == -0.0
true

となるのでこの二つの違いは判別できません。そのため、setrounding をエミュレートできているかをテストする際には、isequal を使いました。

julia> isequal(1.0, 3.0)
false

julia> isequal(1.0, 1.0)
true

julia> isequal(0.0, 0.0)
true

julia> isequal(0.0, -0.0)
false

julia> isequal(-0.0, -0.0)
true

julia> isequal(NaN, NaN)
true

`set_zero_subnormals`

なかなか変更することは少ないとは思われますが、Juliaには非正規化数が0に変換されることを許容することで計算を高速化する(かもしれない)設定が存在します。

Base.Rounding.set_zero_subnormals
Base.Rounding.get_zero_subnormals

上付き丸め、下付き丸めが必要となっている今回のような状況では当然オフになっている必要があります。

RoundingEmulator.jl

add12, mul12

setrounding_raw