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