お話：数値解析 第 9 回
丸め誤差
長田直樹
COBOL などのプログラミング言語、埋め込み式のマイクロプロセッサでは、2 進あるいは 10 進の固数値計算は、無限にある実数を有限個の数 — 手 定小数点数が用いられる。
計算の場合は有限桁の数、計算機による数値計算の1

はじめに場合は有限個の浮動小数点数または固定小数点数 —浮動小数点数3を用いて計算するため、丸め誤差を伴う。
丸め誤差の知識は、数値計算結果の評価に不可欠であり、桁落ちや情報落ちなどの回避にも役立つ。
今回は数値計算における丸め誤差の話をする。

符号と有効数字部 (仮数部ともいう) と指数により表した数を浮動小数点数という。10 進浮動小数点数

6.02 × 1023 の符号は +(正)、有効数字部は 6.02、指数は 23 である。2 進浮動小数点数 1.101(2) × 21 (=3.25(10) ) の有効数字部は 1.101(2) 、指数は 1 である。
計算機では 2 進、10 進、または 16 進の浮動小数
実数を小数で表す際に固定小数点数と浮動小数点 点数が用いられる。固定小数点数に比べ演算に時間数が用いられる。円周率 3.14 は 10 進固定小数点数 がかかるが、表せる数の範囲が広くなるという長所で、アボガドロ定数 6.02 × 1023 は 10 進浮動小数点 があるので、数値計算には通常浮動小数点数が使わ

2固定小数点数数である。計算機で実数を扱うときにも固定小数点 れる。
数と浮動小数点数の２つの方式がある。
整数部と小数部の桁数が決まった表示方法が固定

3.1小数点数である。
整数部 q 桁と小数部 r 桁の 10 進固定小数点数は2 進浮動小数点数

2 進浮動小数点数は
()
b1
bp−1
± b0 + 1 + · · · + p−1 × 2e22

全部で 2 × 10q+r − 1 個 (0 以外は正負 2 つ) ある。
以下では何進数かを明示するときは丸括弧をつけ(1)た下付きの添字を用いる。2 進数は 0 と 1 により表と表される。ここで、bj は 0 または 1、e は整数され、0 と 1 の個数はビットという単位で数える。2で emin 5 e 5 emax である。b0 + b1 2−1 + · · · +進整数と 2 進小数があり

13(10) = 2 + 2 + 2 = 1101(2) ,bp−1 2−p+1 は有効数字部 (仮数部ともいう)、b1 2−1 +· · · + bp−1 2−p+1 は小数部、e は指数という。b0 = 13.25(10) = 21 + 20 + 2−2 = 11.01(2)のとき正規化浮動小数点数 (正規化数と略す) という。

320

ゼロは正規化数ではない。というように表す。1101(2) は 4 ビット 2 進整数、2 進正規化数の個数は 2p (emax − emin + 1) である。

11.01(2) は整数部 2 ビット小数部 2 ビットの 2 進固 最大数は b1 = · · · = bp−1 = 1, e = emax のときで、定小数点数である。
()
1 − 2−p× 2emax = 2emax +1 1 − 2−p
整数部 q ビットと小数部 r ビットの 2 進固定小数−11−2
点数は全部で 21+q+r − 1 個ある。
である。正の最小正規化数は b1 = · · · = bp−1 =0, e = emin のときで 2emin である。
1

3.2

C 言語は歴史的に 0.5 5 有効数字部 < 1 であるので、指数の範囲を IEEE754 より 1 だけ大きくしてIEEE754

PC やワークステーションの 2 進浮動小数点数といる [1]。C の単精度では emin = −125, emax = 128
して IEEE754-1985 が標準的に用いられている。昨で、バイアスは 126 である。たとえば、3.25(10) は年 IEEE754-1985 が 23 年ぶりに改訂された。2008
C の単精度では
年 6 月に IEEE754-2008 が承認され 8 月に出版され()110100た。これまでは、IEEE754-1985 で 2 進浮動小数点++ 3 + 4 + 5 + · · · + 24 × 22
2 222222数、IEEE854-1987 で 10 進浮動小数点数の規格が定められていたが、IEEE754-2008 ではこれらが統合 と表す。ビットパタンは IEEE754 と同一になる。
された。ここでは、IEEE754-1985(IEEE754 と略す)に沿って話していく。

3.3(1) の型の浮動小数点数は符号に 1 ビット (正のとき 0、負のとき 1)、指数部に E ビット、(b0 は 1 と決まっているので省略し) 小数部に b1 , . . . , bp−1 を格IEEE754 での丸め誤差浮動小数点数は有限個であるので、1. 実数を浮動小数点数で近似納する。全体で E + p ビットである。

2. 浮動小数点数同士の四則演算

符
号

指数部

小数部

1Ep−1

(2 つの浮動小数点数を実数と考えたときの和は必ずしも浮動小数点数ではない。)の際に誤差が入る。これらの端数処理に伴う誤差を

指数部は 0 から 2E − 1 まで 2E 通りがあり、(1)E−1の指数 e にバイアス 2− 1 を加えた e + 2

E−1丸め誤差という。

−1数 emin は 2 − 2E−1 、最大指数 emax は 2E−1 − 1 でIEEE754 では、4 種類の丸め方が規定されており、デフォルトは最も近い浮動小数点数に丸める (最も近い値への丸め)。最も近い浮動小数点数が 2 つあるある。したがって、正規化数の指数部は 1 ∼ 2E − 2

ときは、小数部の最下位ビットが 0 になるほうを選を (符号なし)2 進整数で表したものとする。最小指となる。指数部の 0 はゼロと非正規化数、2 − 1 はぶ (IEEE754-2008 では偶数への丸めと呼んでいる)。

無限大と非数に割り当てる。指数部を非負にするこ

実数 x (2emin 5 |x| 5 2emax +1 (1 − 2−p )) を最も近とにより浮動小数点数の大小比較が容易になる。

い値に丸めた浮動小数点数を round(x) と表す。x は∞∑bj 2e , emin 5 e 5 emax
x = ± 1 +2j

EIEEE754 では 4 種類の形式が規定されているが、よく用いられている単精度と倍精度について次表に示す。

j=1

単精度

倍精度

ビット数

3264p2453E811

バイアス1271023

emin−126−1022emax1271023と表せるので、round(x) は、±(1 + b1 2−1 + · · · + bp−1 21−p )2e−1±(1 + b1 2±2e+1+ · · · + bp−1 2

1−p+2(2)1−pe)2(3)(4)のいずれかである。(2) は切り捨て、(3) と (4) は切り上げである。とくに、(4) は繰り上がりが生じている。

たとえば、3.25(10) は単精度では
()
10
100
1 + + 2 + 3 + 4 + · · · + 23 × 21
2 2222

丸め誤差の絶対値の上限は 2e−p である。((1+2−p )2eは 2e に丸められる。このとき丸め誤差は 2e−p である。)

1 より大きい正規化数の最小のものを 1 + eps とす
部 10100000000000000000000(2) である。ビットパタ る。eps をマシンエプシロンという。この定義によンは 0 10000000 10100000000000000000000 となる。 ると IEEE754 では、eps = 21−p となる。
なので、符号 0、指数部 128(10) = 10000000(2) 、小数2

丸めの相対誤差 (= (近似値 − 真値)/真値) の絶対正規化数の加算は次の順序で行う。u, v を正規化たがって、round(x) = x(1 + δ),|δ| 5

正規化数の加算

4

値の上限はマシンエプシロンの半分 12 eps である。し数とする。

1eps
21. u, v の指数部を比較し、u の指数部が小さいととなる。 12 eps = 2−p を丸めの単位という。

きは u, v を入れ替える。

2. v の指数部を u の指数部に揃え、v の有効数字3.4部を変更する。

浮動小数点数の例

分かりやすいように、p = 4, emin = −1, emax = 23. 加算を行う。の場合を考える。すべての正の正規化数を 10 進固4. 正規化を行う。

定小数点数で表したものを次表に示す。有効数字部(SD) は 2 進固定小数点数、指数は 10 進整数で表す。
SD \ 指数−10121.000(2)1.001(2)0.5
0.5625

1.0
1.125

2.0
2.25

4.0
4.5

1.010(2)1.011(2)1.100(2)0.625
0.6875
0.75

1.25
1.375
1.5

2.5
2.75
3.05.05.5
6.0

1.101(2)1.110(2)0.8125
0.875

1.625
1.75

3.25
3.5

6.5
7.0

1.111(2)0.9375

1.875

3.75

7.5

アルゴリズムの詳細は [4, pp.204-207] を見よ。指数部を揃える際に小数部 p − 1 桁では表現できなくなるので、ほとんどの計算機は 2 桁以上の保護桁が用意されている。
例 1 正規化数の加算は ⊕ で表すことにする。3.4 節の正規化浮動小数点数 (丸めは最も近い値への丸め)において (1.001(2) × 21 ) ⊕ (1.001(2) × 20 ) を考える。

+

1.001 × 21
0.1001 × 21

(指数部を揃える)

正規化数は負の数を含め 64 個あり、ゼロに関し対1.1011 × 2であり、1.1011(2) ×21 は両隣の正規化数 1.101(2) ×21

称であるが均等には分布していない。区間 (−0.5, 0.5)と 1.110(2) × 21 の中点なので、小数部の最後のビッには存在せず、[0.5, 1.0] では 0.0625 刻み、[1.0, 2.0]

トが 0 である 1.110(2) × 21 に丸める。したがって、

1

では 0.125 刻み、[2.0, 4.0] では 0.25 刻み、[4.0, 7.5]

(1.001(2) × 21 ) ⊕ (1.001(2) × 20 ) = 1.110(2) × 21

では 0.5 刻みである。正規化数およびゼロの分布を
図 1 に示す。

である。

-2 -1 0123456

保護桁を持たないと1.001 × 217+

図 1: 正規化数とゼロの分布

0.100 × 21(小数第 4 位は落ちる)

1.101 × 21となるので、最も近い値への丸めにならない。

最も近い値への丸めでは、隣り合った正規化数の中点以外は小数部の 2−4 の係数 b4 が 1 のとき切り上げ、0 のときは切り捨てる (零捨一入)。中点のときは b3 が 0 の方に丸める。丸め誤差の絶対値

5

桁落ちと情報落ちの上限は 0.0625(10) × 2 である。マシンエプシロ

5.1

桁落ち

e

ンは 0.125(10) = 1.000(2) × 2−3 で、丸めの単位は近い数同士の引き算をしたとき、有効桁数が少な0.0625(10) である。

くなることを桁落ちという。たとえば、

1.23456 − 1.23460 = −0.000043は有効桁数 6 桁の数同士の引き算であるが、結果は5.2有効桁数が 1 桁になっている。

情報落ち

絶対値の大きな数に絶対値の小さな数を加えると桁落ちは式の変形により避けられることがある。

絶対値の小さな数の情報が反映しない現象を情報落

例 2 |x| が十分小さいときちという。

1.23456 × 103 + 6.54321 × 10−2 を 10 進 6 桁浮動小数点数として計算する。指数を 103 に揃え保護桁を式の通り計算すると桁落ちが生じるが、分子の有 を 2 桁取り有効数字部を計算する。
1.23456
103
理化
x+ 0.0000654 103
f (x) = √1+x+1
1.2346254 103
f (x) =

√1+x−1を行うと桁落ちを避けることができる。
√√
次表に 1 + x − 1 と x/( 1 + x + 1) を示す。有有効数字部の小数第 6 位を四捨五入して1.23456 × 103 + 6.54321 × 10−2 = 1.23463 × 103

効数字部の数字の下付きの添字は同じ数が添字の数だけ続くことを示す。たとえば、4.94 88 = 4.999988が得られる。6.54321 × 10−2 の有効数字部 0.04321である。
√x/( 1 + x + 1)

が結果に反映していない。

1+x−1

1.0 × 10−14

4.8850 × 10−15

4.913 88 × 10−15

例 3 級数

1.1 × 10−145.5511 × 10−15

5.4912 85 × 10−15−14−15√x

1.2 × 105.9952 × 10

5.913 81 × 10−15∑∞−2= π 2 /6 の部分和を
i=1 isν =

ν∑
1i=1

√i2

1 + x − 1 の有効桁数は 2 ∼ 3 桁であるので、 とおき、s8000 を 2 通りの方で計算する。計算は単精√13 ∼ 14 桁桁落ちしている。x/( 1 + x + 1) は 15 桁 度 (10 進 7.2 桁) である。
以上正確なので桁落ちはない。
級数展開

1111(((( 2 + 2 ) + 2 ) · · ·) +)12380002

2x x−28によっても同様に桁落ちを避けることができる。
f (x) ∼により計算すると 1.6447253 となる。ν = 4096 のとき sν = 1.6447253 である。これは、桁落ちは |x| が 0 に近いとき、

1 − cos x,

11= 2−24 = eps,40962
2sin(a + x) − sin a

1 < sν < 2であるので、ν = 4097 のとき 1/ν 2 は結果に反映していない。一方、加える順序を逆にしてなどでも生じる。それぞれ、

1 − cos x = 2 sin2 x2((((sin(a + x) − sin a = 2 sin x2 cos(a + x2 )1111+)+) · · ·) + 2 )
80002
79992
79982
1により計算すると 1.6448090 と情報落ちは回避できる。と変形することにより桁落ちは避けられる。桁落ち防止の変形の多くは、不定形の極限を求めるときに級数の和を高精度に計算する方法にカハンの加算用いる技法と共通している。

アルゴリズムがある。浮動小数点数 x1 , · · · , xν の和x → ±∞ のときにも桁落ちが生じることがある。 ∑νi=1 xi を次のアルゴリズムで計算する。[2, pp.83ex
1288]
1 + tanh x =, 1 − tanh x =cosh x
cosh2 xを用いた桁落ち防止の変形を、二重指数関数型数値積分公式 (2008 年 7 月号) で行った。前者は x → −∞のとき、後者は x → ±∞ のとき桁落ちが生じる。

4s=0s = t + y の際に情報落ちとなる x(2) + e(2) は e に格納される。y = xi + e における情報落ちのため、ee=0
for i = 1 : νtmp = sは正確な補正ではない。
筆者の同僚にベクトルの内積の高精度計算 [5] な

y = xi + es = tmp + yどの研究をしている荻田武史さんがいる。彼に尋ねた所、カハンのアルゴリズムにより倍精度 (64 ビッtmp = tmp − se = tmp + yendト) で計算する際、浮動小数点レジスタが 64 ビットの場合はうまく働くが、浮動小数点演算の誤差を利用しているため浮動小数点レジスタが 80 ビットの場合は問題が生じることがあるという。
ためしに、倍精度で定理 1 (クヌース [4]) カハンの加算アルゴリズムで計算した和 s は

s=

ν
∑

1111
(((( 4 + 4 ) + 4 ) · · ·) +)123160004|δi | 5 2u + O(νu2 )(1 + δi )xi ,をカハンのアルゴリズムにより計算してみた。浮動i=1小数点レジスタが 80 ビットである Vine Linux 3.3を満たす。ここで、u = 2−p は丸めの単位である。

上 gcc 3.3.2 と Cygwin 上 gcc 3.4.4 を用いたが、いずれもカハンのアルゴリズムは正常に働いた。

¤証明 [4, p.230,p.537][2, p.85] を見よ。
例 4 カハンの加算アルゴリズムにより61111s8000 = (((( 2 + 2 ) + 2 ) · · ·) +)12380002関孝和の円周率の計算の誤差
連載の第 5 回 (2008 年 9 月号) で関孝和の円周率を単精度で計算すると s8000 = 1.6448091 となり、情の計算を取り上げた。直径 1 の円に内接する正 2ν 角

報落ちは回避できる。

形の周の長さを sν とする。
「s15 , s16 , s17 の小数第 17

カハンの加算アルゴリズムにおいて、s =位以下は桁落ち (近い数同士の引き算により有効桁∑i−1j=1 xj

数が少なくなる現象) などのため正しくない。」と書に x = xi を加えるところの有効数字部の変化を図いたが、厳密にいうと桁落ちとも情報落ちとも異な解する。有効数字部を p ビットとし、s と x の指数り、固定小数点数計算で小さい数を扱ったことによの差が q であるとする。図 2 では指数部を揃え、有る有効桁数の減少が原因である。

効数字部に 2p ビット使っている。s の上位の q ビッ
トを s(1) とし、下位の p − q ビットを s(2) とする。

x の上位の p − q ビットを x(1) 、下位の q ビットをx(2) とする。e の上位の p − q ビットを e(1) 、下位のq ビットを e(2) とする。桁数を揃えるため、必要なときは上位や下位ビットを 0 で埋める。

ビット数

qp−q

t=s(1)(2)sy = xi + e
s=t+yt=t−ss+ e(1)

(2)

(1)

s−x+y(1)e=t+yp−qs(1)x(1)qx(2) + e(2)

図 3:『括要算法』(東北大学・岡本刊 089)

−e(1)

0(2)

x+ e(2)

関の没後に出版された『括要算法』[7] に計算結果0(図 3) が記載されている。左側のページの 3 行目までは周の計算、「第二 求定周」ではエイトケン ∆2 法

図 2: カハンの加算アルゴリズム

5参考文献と呼ばれている加速法により円周の長さを定めている。一絲 (し) は一尺 (約 30.3cm) の 105 分の一、一塵 (じん) は一尺の 1010 分の一である。直角三角形[1] 日本工業標準調査会、プログラム言語 C JIS Xの斜辺を弦、直角を挟む短い辺を勾 (こう)、長い辺3010:2003(”ISO/IEC 9899:1999 - Programming Language C”の日本語訳。)

を股 (こ) という。

0.90.8

http://www.jisc.go.jp/C1AD[2] N.J. Higham, Accuracy and Stability of Nu-B

merical Algorithms, 2nd. ed., SIAM, 2002.

0.70.6[3] IEEE, Draft Standard for Floating-Point
Arithmetic P754, 2006.

0.50.4

http://isavhdlimp.sourceforge.jp/pdf/IEEE754r-2006-10-04.pdf

0.30.2

O0.1[4] D.K. クヌース著、有沢・和田監訳、The Art ofComputer Programming, Vol.2, Seminumerical Algorithms, 3rd. ed. 日本語版、ASCII, 2004

0−0.6

−0.4

−0.2

00.20.4

図 4: 正 2xν 角形の勾股弦0.6

図 4 において、AB は正 2ν−1 角形の弦、CD、AD とAC はそれぞれ正 2ν 角形の勾、股、弦である。
規則性を見て取れるように、ν = 17(一十三萬一千[5] 荻田武史、品質を落とさない数値計算法、数学セミナー、2008 年 11 月号、15-19

零七十二角) の場合に一尺を 1 とする固定小数点数[6] 長田直樹、関孝和の円周率の計算、京都大学数理解析研究所講究録、近刊
http://www.cis.twcu.ac.jp/~osada

で表した値を次表に示す。
固定小数点数表示
勾0.00000 00005 74486 5862 強/rims2008.pdf

股弦

0.00002 39684 49801 5334 強
0.00002 39684 49808 4182 強[7] 関孝和、括要算法、東北大学和算ポータル、岡周3.14159 26532 88992 7759 弱

本刊、請求番号 089

http://www2.library.tohoku.ac.jp/wasan/
関は小数第 20 位を丸めて第 19 位まで (整数部 1(おさだ なおき/東京女子大学)

桁、小数部 19 桁の固定小数点数により) 表している
ことが分かる。
小数第 d 位まで計算したとすると、勾の有効桁数は d − 9 である。関は√22弦= 勾 +股 ,周 = 2ν 弦によって周を求めたので、周の有効桁数も d − 9 程度となる。周の有効桁数は 17 であるので、勾の有効桁数も 17 桁前後と考えられる。つまり関の誤差は、有限桁の固定小数点数計算において計算の途中で絶対値の小さい数が現れたことによる有効桁数の減少が原因、ということになる。
有効桁数の減少を手がかりにして、関が図 3 の値を小数第 26(= 17 + 9) 位までの計算により求めたことが示されている [6]。

6