第7回：■ 浮動小数点数

■ 浮動小数点数

▼ 正規化数、副正規化数

浮動小数とは、$0.12$ の代わりに $1.2 \times 10^{-1}$ のように表示することである。

10進数の浮動小数は

\[\pm{\left(d_0.d_1d_2\cdots \right)}_{10}\times 10^{e}\]

のように表される。$\times$の前までの${\left(d_0.d_1d_2\cdots \right)}_{10}$ の部分は仮数部と呼ばれる。添字の$10$は10進数を意味し、$d_0, d_{1}, \cdots $ は $0,1,\ldots,9$ までの数字である。$\times$ の後ろの $10^{e}$ は指数部と呼ばれる。

2進数の浮動小数は

\[\pm{\left(b_0.b_1b_2\cdots \right)}_{2}\times 2^{e}\]

のように表される。ここで、$\times$の前までの ${\left(b_0.b_1b_2\cdots \right)}_{2}$ の部分は仮数部と呼ばれる。添字の$2$は2進数を意味し、$b_0, b_{1}, \cdots $ は $0$ または $1$ の数字である。$\times$ の後ろの $2^{e}$ は指数部と呼ばれる。

bit (binary digit)とは、2進数の一桁のことである。

本文で用いる浮動小数点数は Float64 型である。

julia> typeof(1.0)
Float64

Note

以下で、浮動小数点数の2進数による表現を詳しく説明するが、丸暗記する内容ではない。しかし、計算機内部の小数が「有限桁」で行われることは、計算機による数値計算では常に意識すべきである。

Float64 型は、「IEEE754標準倍精度浮動小数点数」に基づき、符号部 1 bit、指数部 11 bit 仮数部 53 bit から構成される。ただし、以下のように先頭の 1 bitを固定し、仮数部の 52 bit のみをデータとして採用するため、 2進数の並びは 1+11+52 = 64 bit である。

Float64は、正規化数、副正規化数、数でない数の３種類からなりたっている。

正規化数は、$b_{0} = 1$として、 $\pm{\left(1.b_1b_2\cdots b_{52}\right)}_{2}\times 2^{e}$ のように表すものである。ただし、指数は $−1022 \le e \le 1023$ の範囲である。仮数 ${\left(1.b_1b_2\cdots b_{52}\right)}_2$ は 1以上で2を超えない範囲の小数となる。

正規化数で表すことができない、絶対値が小さい浮動小数は副正規化数で表わされる。

副正規化数は、$b_{0} = 0$, $e=−1023$ として、 $\pm{\left(0.b_1b_2\cdots b_{52}\right)}_{2}\times 2^{e}$ のように表すものである。仮数部 ${\left(1.b_1b_2\cdots b_{52}\right)}_{2}$ は 0以上で1を超えない範囲の小数となる。

「数でない数」は、■ 0による除算で、既に説明した。 Inf, -Inf, NaN の３つである。

Float64で表すことができる、絶対値が最も大きい数は、正規化数の $2^{1024}≃1.798\times10^{308}$ である。絶対値が最も小さい数は、副正規化数の $2^{−1022}≃2.225\times10^{−308}$ である。

これらは、関数 realmax, realminでそれぞれ得られる。

julia> realmax(Float64)
1.7976931348623157e308

julia> realmin(Float64)
2.2250738585072014e-308

丸め

小数 $0.2$ は $0.2 = \frac{1}{5} = \frac{1}{{101}_{2}}$ となるが、$1$ を ${101}_{2}$ で割り切ることはできない。$0.2$ を2進数で表すと

\[{0.00110011001100\cdots}_{2}\]

のようになる。すなわち、$1100$ の並びが無限に続く循環小数となる。

また、小数 $0.1$ は $0.1 = \frac{1}{5\times 2} = \frac{1}{{101}_{2}} \times 2^{-1}$ であるから、$0.1$ を2進数で表すと（上を1桁ずらして）

\[{0.000110011001100\cdots}_{2}\]

のようになる。これも、$1100$ の並びが無限に続く循環小数となる。

「有限桁の小数」で表すことができない「循環小数」を、 Float64型で表現するとき、その仮数の末尾に近いの桁を修正する操作を行う場合がある。この操作を「丸める」という。

「丸め」られた浮動小数の計算は、筆算とは違う結果となる場合がある。例えば、

julia> 0.1+0.2
0.30000000000000004

julia> 0.1+0.2 == 0.3
false

筆算の結果は $0.3$であるが、計算結果は 0.30000000000000004 と異なってしまう。

別の例として、$0.1$を 10回足した結果は

julia> s=0
0

julia> for i in 1:10
         s += 0.1
       end

julia> @show s
s = 0.9999999999999999
0.9999999999999999

julia> s == 1.0
false

0.9999999999999999 となり、$1.0$ にはならない。

このような、「丸め」を原因とする、正しい値からの「ずれ」を「丸め誤差」と呼んでいる。

▶︎ 小数を2進数へ変換する

\[x=f_{1}2^{-1} + f_{2}2^{-2} + \cdots\]

(正の)小数を2進数に変換するには、小数を2倍しその整数部分を取り出すことを、繰り返し行えばよい。

小数 0.2を、2進数で表示すると循環小数になる。 1100 のパターンが繰り返し現れる。

julia> x=0.2
0.2

julia> for i=1:50
           q=floor(x/2)
           print(Int64(q))
           x -= q*2
           x *= 2
       end
00001100110011001100110011001100110011001100110011

上の結果の最初の桁は、$2^1$ の桁に相当する。すなわち、小数点は、２つ目の数字の後ろに位置する。

◀︎ 練習：有限小数・循環小数

0.5以下の正の小数をいくつかを選び、これらを2進数に直してみよ。有限小数か循環小数かを判定せよ。

例: 0.1, 0.2, 0.25, 0.3, 0.5

さらに、5つ程度の例を加えよ。

■ 加減算における桁落ちと情報落ち

加算と減算は、小数点の位置を合わせて計算されるが、桁数が有限であることから、正しい得られない場合がある。その原因のうち「桁落ち」と「情報落ち」の二つの現象が知られている。

■ 桁落ち

「桁落ち」は、互いに非常に近い二つの数 $x, y$ に対して、減算 $x-y$を行うと、結果の有効桁数が大きく減少する現象である。

例えば、有効桁数が4桁の二つの数の引き算の例を見よう。

julia> 2.345 - 1.233
1.112

julia> 1.234 - 1.232
0.0020000000000000018

前者の結果は 4桁の有効桁数を保っているのに対して、後者の結果は 1桁の有効桁数になってしまう。（末尾の 18 は丸め誤差である)。

式を変形して、互いに近い数同士を引くことを回避できる場合がある。下の例を参考にせよ。 → ▶︎ 2次方程式

■ 情報落ち

「情報落ち」は、絶対値が大きく異なる数を加減算すると、小さい桁の精度が失われる現象である。

例えば、3つの数 $x = 14\times 10^{-17}$, $y = 24\times 10^{-16}$, $z = 1$ を、この順番で加えた結果と、逆の順番で加えた結果を比較しよう。

julia> x=14e-17
1.4e-16

julia> y=24e-16
2.4e-15

julia> z=1
1

julia> xyz=(x+y)+z
1.0000000000000024

julia> zyx=(z+y)+x
1.0000000000000027

筆算による正しい値は 1.00000000000000253 であるが、後者の和よりも前者の和が、正しい値に近い。

後者の和が誤差を大きく含んだのは、和 $z+y$ の段階で、有効桁数をほぼ使い切ったからである。

julia> zy=z+y
1.0000000000000024

julia> nextfloat(zy) # z+y の「隣りの」正しく表される数
1.0000000000000027

一般に、大きさの異なる数同士を加減算する場合には、絶対値が小さいものから計算を進めたほうがよい。

ここで見たように、有限桁数の浮動小数点数の加減算は「結合則」を満たさない。

\[(x+y)+z \neq x+(y+x)\]

■ 等比級数

Base.logspace — Function

関数 logspace(a,b,n=50) は、等比級数を作る方法の一つである。初項 $10^{a}$ から始めて $10^{b}$ で終わる等比級数 (要素は$n$個)となる ■ ベクトルを作る。

以下は、初項 $10^{0}=1$ から始めて、$10^{3}=1000$ で終わる等比級数 (要素は $4$ 個 ) を作る。すなわち、$1, 10, 100, 1000$ が作られる。

julia> logspace(0,3,4)
4-element Array{Float64,1}:
    1.0
   10.0
  100.0
 1000.0

以下は、初項 $10^{-3}$ から始めて、$10^{0}=1$ で終わる等比級数 (要素は $4$ 個 ) を作る。すなわち、$0.001, 0.01, 0.1, 1.0$ が作られる。

julia> logspace(-3,0,4)
4-element Array{Float64,1}:
 0.001
 0.01
 0.1
 1.0

Note

logspace = logarithmically spaced numbers

2次方程式 $x^2-bx+c=0$ の解は、解の公式から、判別式 $d=b^2-4c$を用いて、 $\begin{align}x_1&=\frac{b+\sqrt{d}}{2}=\frac{b+\sqrt{b^2-4c}}{2}\\ x_2&=\frac{b-\sqrt{d}}{2}=\frac{b-\sqrt{b^2-4c}}{2}\end{align}$ であるが、$b$ と $\sqrt{d}$が同程度のとき $x_2$ は「桁落ち」しやすい。

そこで、$(b+\sqrt{b^2-4c})$ を分母分子に掛けて

\[x_{21} = \frac{2c}{b+\sqrt{b^2-4c}}=\frac{c}{x_1}\]

のように変形してから計算する。最後の項は、解と係数の関係 $x_1x_2=c$ である。

\[\begin{align*} (x-x_1)(x-x_2)&=x^2-(x_1+x_2)x+x_1x_2 = x^2-bx+c,\\ b&=x_1+x_2, \\ c&=x_1x_2 \end{align*}\]

▶︎ 2次方程式：計算の例

実例で見てみよう。

小さい正の数 hを用いて、$\alpha = 100+h$ と $\beta = 1+h$ を解とする2次方程式を作る。解と係数の関係から、上の方程式において $b = \alpha + \beta$, $c=\alpha\beta$ と定めればよい。

h=logspace(-12,-1);
alpha=100+h
beta=1+h;
c=alpha .* beta;
b=alpha .+ beta;

解の公式から、「大きい方の解」 x1を計算する。 x2sは解の公式から求めた「小さい方の解」、である x2vは解と係数の関係から求めた「小さい方の解」

d=b.*b-4c;
x1=(b+sqrt.(d))/2;
x2s=(b-sqrt.(d))/2;
x2v=c./x1;

「大きい方の解」について、正しい解との差をプロットしてみる。

using PyPlot
plot(h, x1-alpha, ".")
xlabel("h")
ylabel("x1-alpha")
xscale("log")

「小さい方の解」について、正しい解との差をプロットしてみる。

plot(h, x2s-beta,".",label="x2s")
plot(h, x2v-beta, "o",label="x2v")
xlabel("h")
ylabel("x2-beta")
xscale("log")
legend()

「小さい方の解」について、正しい解との差の絶対値(残差)をプロットしてみる。

plot(h, abs.(x2s-beta),".",label="x2s")
plot(h, abs.(x2v-beta), "o",label="x2v")
xlabel("h")
ylabel("abs(x2-beta)")
xscale("log")
ylim(1e-18,1e-13)
yscale("log")
legend()

解の公式から求めた「小さい方の解」の残差が「あばれる」のに対して、解と係数の関係から求めた小さい方の解」の残差が「一定」である様子が見れる。

▶︎ 数値微分

\[\frac{df(x_0)}{dx} = \lim_{h \longrightarrow 0}\frac{f(x_0+h)-f(x_0)}{h}\]

関数 $y=x^2$ の $x=1$ における微分係数を、上の定義により求めよう。求まるべき値は $2$ であるが、$h$ を小さくすると $2$ の上下に暴れてしまう。

using PyPlot
h=logspace(-18,-8,100)
d=( (1+h).^2 - 1) ./ h
plot(h,d, ".")
ylim(5e-1,3e0)
yscale("log")
xscale("log")

今度は、関数 $y=x^n$, ($n=1,2,3$) の $x=1$における微分係数を、上の定義により求めよう。求まるべき値は $n$ であるが、$h$ を小さくすると $n$ の上下に暴れてしまう。

using PyPlot
h=logspace(-18,-8,100)
for n=1:3
    d=( (1+h).^n - 1) ./ h
    plot(h,d, ".", label="y=x^"*string(n))
end
xlabel("h")
ylabel("d")
yscale("log")
xscale("log")
legend()

以上の誤差も、非常に近い二つの数字を減じたときに現れる「桁落ち」の現象である。 ▶︎ 2次方程式とは異なり、うまく回避する手段はない。$h$を小さく取りすぎないように注意する。

▶︎ 練習・数値微分

以下の関数の、指定された座標での微分係数を、上の例と同様に求めてみよ。

指数関数 $y = \exp{x}, \; x = 0$
対数関数 $y = \log{x}, \; x = 1$
対数関数 $y = \log\left(1+x\right), \; x = 0 \;$ ▶ 注: 関数 Base.log1p — Function を用いよ。
三角関数 $y = \sin{x}, \; x = 1$。▶ 正しい微分係数は 0.540302305868140 である。

■ 近似比較演算子 `isapprox`

条件式 x == 1 は、数 x が 1 と完全に一致することを判定するので、数 x が丸め誤差を含むような場合に用いるのに適さない。

その代わりに、丸め誤差の基準を適当な数、例えば、$10^{-6}$ をとって、条件式 abs(x-1) < 1e-6 をもって、数 $x$が $1$ に非常に近いことを判定するのが常套手段である。

julia> x=1+1e-8
1.00000001

julia> x == 0
false

julia> abs(x-1) < 1e-6
true

Base.isapprox - Function

Julia には、数 a と b がほぼ等しいことを判定する近似比較演算子 isapprox(a,b) が用意されているので、必要に応じて用いるとよい。a と b との丸め誤差の程度を考慮して、比較を行う便利な関数である。

julia> 0.1+0.2 == 0.3
false

julia> isapprox(0.1+0.2, 0.3)
true

■ 数でない数の判定

Numeric Comparisons

■ 0による除算で社迂回したように、 IEEE754規格の浮動小数点数は、「数でない数」NaN, Inf, -Inf の3つを含んでいる。これらを判定する関数が用意されている。

julia> for x in [0,1,Inf,NaN,NaN]
         println()
         @show isfinite(x)
         @show isinf(x)
         @show isnan(x)
       end

isfinite(x) = true
isinf(x) = false
isnan(x) = false

isfinite(x) = true
isinf(x) = false
isnan(x) = false

isfinite(x) = false
isinf(x) = true
isnan(x) = false

isfinite(x) = false
isinf(x) = false
isnan(x) = true

isfinite(x) = false
isinf(x) = false
isnan(x) = true

★今回のまとめ

浮動小数点数
有限小数・循環小数
加減算における桁落ち・情報落ち
近似比較演算子
等差数列・等比数列
数値微分
数でない数