第一回課題解説

大石進一

2001 年 5 月 4 日

全体によくできていた。採点がとても楽しく感じられた。以下，問題を示し，各問ごとに解説を加えることにしよう。

1 問題

1. 教科書「数値計算」 (p.11 問題 1.) 1 から 4 を解け。以下にその問題を示す。

1-1. 倍精度浮動小数点数

▽(1/10)

= +0.11001100110011001100110011001100110011001100110011001_{2} × 2^{−11_{2}}
= 0.09999999999999999167332731531132594682276248931884765625

を 10 進数の分数で表せ。 1-2. 使っている OS の下での C 言語の printf がどんな中身になっているか調べよ。 1-3. |x| が 0 に近いとき、次の式を不動小数点数を使って評価する:

\frac{1}{1 + 2x} − \frac{1 − x}{1 + x}

この式と等価な式で桁落ちを少なくする式としてどんなものがよいか考えよ。 1-4. C 言語において、 10 進数から 2 進数に変換する際にコンパイラの最適化の影響を避け

るためには、 9 ページのプログラムにおいて、 double という宣言を volatile double に置き換えればよい。 volatile の仕様についてマニュアルを調べて、その理由を述べよ。

2. IEEE754 の規格についてまとめよ。

3. Intel 上の C での IEEE754 の double の 64bit を書き出すプログラム [1] がある。これ、または、自作でプログラム作り、実際に動かして、 Intel の現在の CPU の nearest mode に IEEE754 の規格違反があるかどうか論ぜよ。

2 解答 1-1

易しい問題で，幾つかの解答パターンがあった。その幾つかを示そう。解答パターン 1 (河邊昌彦君他)

問題の数は、

+0.11001100110011001100110011001100110011001100110011001_{2} × 2^{−11_{2}}
= 11001100110011001100110011001100110011001100110011001_{2} × 2^{−3} × 2^{−53}
= 11001100110011001100110011001100110011001100110011001_{2} × 2^{−56}
= \frac{11001100110011001100110011001100110011001100110011001_{2}}{2^{56}}

と変形できる。この分数の分子は 2 を因数として持たないので、全体としては規約分数である。従って、 10 進数の分数で表すと、

\frac{7205759403792793}{72057594037927936}

となる。解答パターン 2 (川岸寛之君他) 1/10 は 0.1 = 0.000\overline{1100}_{2} = 0.\overline{1100}_{2} × 2^{−11_{2}} で表される。倍精度の仮数部が 53bit であること

より、

仮数部 = 0.1100_{2} + 0.1100_{2} × 0.0001_{2} + ⋯ + 0.1100_{2} × 0. \overbrace{00 ⋯ 00}^{48−1=47} 1_{2} + 0. \overbrace{00 ⋯ 00}^{52} 1_{2}
= (\frac{1}{2} + \frac{1}{4}) + (\frac{1}{2} + \frac{1}{4}) × \frac{1}{16} + ⋯ + (\frac{1}{2} + \frac{1}{4}) × (\frac{1}{16})^{12} + 2^{−53}

= \frac{\frac{3}{4} × (1 − (\frac{1}{16})^{13})}{1 − \frac{1}{16}} + 2^{−53}
= \frac{7205759403792793}{9007199254740992}
指数部 = 2^{−11_{2}} = 2^{−3} = \frac{1}{8}  (1)

よって、 1/10 = 仮数部 × 指数部 = 7205759403792793/72057594037927936

3 解答 1-2

この問題はかなり難しく感じられたようだ。まず，題意がよく把握できなかったようだ。実は， printf を改良して，丸めを上，下，最近点としたときに，指定した丸めで double の 10 進数表示をする改良された printf を作れという問題の第一段階として，この問題が出されている。

次は河邊昌彦君の解答であり，非常に優れている。

3.1 概要

C 言語の処理系の中で、著名でかつソースがフリーである、 GCC[3] の printf を調べる。使用した GCC のバージョンは 2.95.3 であり、アーカイブファイルを解凍すると、 libio 以下に入出力関係の C ランタイムライブラリのソースがある。

printf は多機能な関数であり、ソースも必ずしも綺麗であるとは言えない。従って、浮動小数点数に関する部分、特に 2 進から 10 進への変換に絞って調べていく。

3.2 printf の書式指定

printf の書式指定を浮動小数点数に限ってまとめる。

%[flags] [width] [precision] type

各フィールドの意味は以下の通りである。 flags

フラグ意味
-  指定されたフィールド幅に結果を左詰めする。
+  出力値の前に+または-の符号を付ける。
0  最小幅まで 0 が付加される。
(空白) 出力値の前に空白が付く。
#  e, E, f の各書式と一緒に使うと、必ず出力値に強制的に小数点が入る。

g または G の書式と一緒に使うと、必ず出力値に強制的に小数点が入り、後続する 0 が切り捨てられない。

width

表示する最小桁数を指定する。 precision

精度を指定する。 type によって意味が変わる。

type 意味
e, E, f 小数点以下の表示桁数を指定する。
g, G 表示する最大有効桁数を指定する。

type パラメータのデータ型は全て double である。

文字出力形式
e  [-]d.dddd e sign ddd
E  [-]d.dddd E sign ddd
f  [-]dddd.dddd
g  e 又は f のうち、短くなる方の形式。
G  E 又は f のうち、短くなる方の形式。

3.3 ファイル構成

printf に関係するファイルをまとめておく。全てのファイルは libio の中にある。

ファイル名行数内容
iostdio.h  83 printf の定義。 _IO_printf に define されている。
ioprintf.c  33 _IO_printf の定義。 _IO_vfprintf が呼ばれている。
iovfprintf.c 185 _IO_vfprintf の定義。

                 浮動小数点数に対しては_IO_outfloat 又は__cvt_double が呼ばれる。
                 (_IO_USE_DTOA の定義に依る)

            701 __cvt_double の定義。

outfloat.c  36 _IO_outfloat の定義。 _IO_dtoa を呼んでいる。
floatconv.c 1753 _IO_dtoa の定義。

libio.h  87 _IO_USE_DTOA の定義。 _G_HAVE_PRINTF_FP の値に依る。

gen-params 720 _G_HAVE_PRINTF_FP の定義の生成。

3.4 実装

3.4.1 printf

iostdio.h(83)

[figure]

#define printf _IO_printf

[figure]

このように、 printf は_IO_printf に define されている。 stdio/stdio.h には printf の外部参照宣言があるが、恐らくそれ以前に iostdio.h がインクルードされているのであろう。

3.4.2 IO printf

ioprintf.h(33)

[figure]

int
_IO_printf
#ifdef __STDC__
  (const char* format, ...)
#else
(format, va_alist) char *format; va_dcl
#endif
{
  int ret;
  va_list args;
  _IO_va_start(args, format);
  ret = _IO_vfprintf(_IO_stdout, format, args);
  va_end(args);
  return ret;
}

[figure]

_IO_printf ではパラメータの調節だけを行って、主な処理は_IO_vfprintf で行っている。実際に、 printf 系の関数は全て_IO_vfprintf で主な処理を行っているようである。

3.4.3 IO USE DTOA

libio.h(86)

[figure]

#if !_G_HAVE_PRINTF_FP
# define _IO_USE_DTOA 1
#endif

[figure]

_IO_USE_DTOA は、 _IO_vfprintf の中で、浮動小数点数の出力にどの関数を使うかに影響を与

える。 _G_HAVE_PRINTF_FP が 0 に定義されていれば、 _IO_USE_DTOA が定義されることになる。

3.4.4 G HAVE PRINTF FP

gen-params(712)

[figure]

# A little test program to check if __printf_fp is available.
cat >dummy.c <<EOF
int main()
{
    return __printf_fp ();
}
EOF
if ${CC} dummy.c >/dev/null 2>&1 ; then
  echo "#define ${macro_prefix}HAVE_PRINTF_FP 1"
  echo "#define ${macro_prefix}HAVE_LONG_DOUBLE_IO 1"
else
  echo "#define ${macro_prefix}HAVE_PRINTF_FP 0"
  echo "#define ${macro_prefix}HAVE_LONG_DOUBLE_IO 0"
fi

[figure]

このシェルスクリプトによって、 __printf_fp という関数があるかどうかがテストされている。その結果、コンパイルに成功したならば、 _G_HAVE_PRINTF_FP を 1 に定義するプリプロセッサ文が生成され、失敗したならば 0 に定義するプリプロセッサ文が生成される。 GCC を使って GCC をコンパイルする場合は、このテストに失敗すると思われるので、 _G_HAVE_PRINTF_FP は 0 に定義されることになる。 3.4.5 IO vfprintf

iovfprintf.h(185)

[figure]

  int
  _IO_vfprintf (fp, fmt0, ap)
       _IO_FILE *fp;
       char const *fmt0;
       _IO_va_list ap;
  {

[figure]

_IO_vfprintf の定義がここから始まっている。処理の流れは、フォーマット文字列を走査して、書式指定ならパラメータを書式指定に沿って出力、それ以外ならフォーマット文字列をそのまま出力、となっている。

浮動小数点数の出力は、 _IO_USE_DTOA が定義されているかどうかで、 _IO_outfloat を使う処理と、 __cvt_double を使う処理に分かれる。前述のように、 _IO_USE_DTOA は定義されているので、 _IO_outfloat の方が使われることになる。

iovfprintf.h(363)

[figure]

case 'e':
case 'E':
case 'f':
case 'F':
case 'g':
case 'G':
_double = va_arg(ap, double);
#ifdef _IO_USE_DTOA
{
    int fmt_flags = 0;
    int fill = ' ';
    if (flags & ALT)
        fmt_flags |= _IO_SHOWPOINT;
    if (flags & LADJUST)
        fmt_flags |= _IO_LEFT;
    else if (flags & ZEROPAD)
        fmt_flags |= _IO_INTERNAL, fill = '0';
    n = _IO_outfloat(_double, fp, ch, width,
     prec < 0 ? DEFPREC : prec,
     fmt_flags, sign, fill);
    if (n < 0)
        goto error;
    ret += n;
}
/* CHECK ERROR! */
continue;
#else

[figure]

_double には書式指定に対応するパラメータが入る。 fp は_IO_vfprintf に渡されたファイルポインタである。 ch は書式指定の型フィールド文字 (“e” 等) である。 width は文字幅指定であり、最小文字列数の指定である。 prec は精度指定であり、何も指定されていないか、或いは負の値が指定された時は-1 である。

従って、そのような場合は DEFPREC(6 に define されている) が精度として使われる。

ALT というフラグは、書式指定のフラグにおいて “#” が指定されていた時にセットされる。同様に、 LADJUST というフラグは “-” が指定されていた時に、 ZEROPAD は “0” が指定されていた時に、セットされる。

sign には符号の付け方を表す文字が格納されており、デフォルトで “\0”、 “ ” を指定すると “ ”、 “+” を指定すると “+” になる。

また、 _IO_outfloat は成功すると書き込んだ文字数を、失敗すると EOF を返す。 EOF は (-1) に define されているので、 _IO_outfloat に失敗すると、エラー処理にジャンプすることになる。 3.4.6 IO outfloat

outfloat.c(36)

[figure]

int
_IO_outfloat (value, sb, type, width, precision, flags, sign_mode, fill)
double value;
_IO_FILE *sb;
int type;
int width;
int precision;
int flags;
int sign_mode;
int fill;
{

[figure]

_IO_outfloat の定義がここから始まっている。渡されるパラメータは、前述の通りである。

outfloat.c(71)

[figure]

switch (type)
{
case 'f':
mode = 3;
break;
case 'e':
case 'E':
exp = type;
mode = 2;
if (precision != 999)
  precision++;  /* Add one to include digit before decimal point. */
break;
case 'g':
case 'G':
exp = type == 'g' ? 'e' : 'E';
if (precision == 0) precision = 1;
if (!(flags & _IO_SHOWPOINT))
  skip_zeroes = 1;
type = 'g';
mode = 2;
break;
}

[figure]

パラメータのセットアップが行われている。

f e E g G
mode 3 2 2 2 2
type ’f’ ’e’ ’E’ ’g’ ’g’
exp  ’e’ ’E’ ’e’ ’E’

ただし、 _IO_vfprintf から呼ばれた時は、 precision は 999 になることはない。

outfloat.c(96)

[figure]

    p = _IO_dtoa(value, mode, precision, &decpt, &sign, &end);

[figure]

_IO_dtoa を呼んで、浮動小数点数を文字列に変換している。 _IO_dtoa のコメントの一部を引用すると、 mode の解釈は以下の通りである。

     2 ==> max(1,ndigits) significant digits.  This gives a
             return value similar to that of ecvt, except
             that trailing zeros are suppressed.
     3 ==> through ndigits past the decimal point.  This
             gives a return value similar to that from fcvt,
             except that trailing zeros are suppressed, and
             ndigits can be negative.

そして、 p には末尾の 0 を除いた有効数字の文字列へのポインタが格納される。また、 decpt, sign には ecvt, fcvt と同様に、小数点と符号の位置が、 end には p の指す文字列の末尾へのポインタが返ってくる。ただし、 value が ±∞, NaN の場合は decpt==9999 となる。

これ以降の処理は、変換の結果からフラグやフォーマットに沿って、符号、仮数部、指数部、小数点、パディングを表示している。

3.4.7 IO dtoa

floatconv.c(1753)

[figure]

   char *
  _IO_dtoa
  #ifdef KR_headers
          (d, mode, ndigits, decpt, sign, rve)
          double d; int mode, ndigits, *decpt, *sign; char **rve;
  #else
          (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
  #endif
  {

[figure]

この関数の定義はかなり長いので、行数のみで引用することにする。

±∞, NaN, 0 の場合 (1833–1868)

まず、特殊な場合に関しての処理を済ませている。これらの場合は、適切な文字列 (“Infinity”, “NaN”, “0”) を返すだけである。なお、符号の処理はこれらの場合にも行われている。 10 を底とした指数部の計算 (1870–1937)

浮動小数点数の 2 を底としている指数部を、 10 を低とする数に変換している部分である。即ち、元の浮動小数点数を d = m × 2^{e} とすると、

log_{10} d = log_{10}(m × 2^{e})
= log_{10} m + e log_{10} 2
= \frac{1}{log 10}(log m + e log 2)

である。ここで、 1 ≤ m < 2 であるので、 log m を 1.5 まわりにテイラー級数展開すると、

log m = log 1.5 + \frac{m − 1.5}{1.5} + o(m − 1.5)

となる。従って、

log_{10} d = \frac{1}{log 10} (log 1.5 + \frac{m − 1.5}{1.5} + o(m − 1.5) + e log 2)
= \frac{log 1.5}{log 10} + \frac{m − 1.5}{1.5 log 10} + e\frac{log 2}{log 10} + \frac{o(m − 1.5)}{log 10}
≈ 0.1760912590558 + (m − 1.5) × 0.289529654602168 + e × 0.301029995663981

である。そして、この右辺を

0 ≤ ∃m′ < 10
: m × 2^{e} = m′ × 10^{e′}

となるような整数 e′ に近似する。

10 進への変換 (2054–2093)

m′ は 0 ≤ m′ < 10 を満たすので、その整数部は高々1 桁である。従って、 m′ の整数部を取り除きながら、 10 を掛けていけば各桁の数 m′_{i} を得ることができる。そして、 m′_{i} と e′ が分かれば、 10 進数として表すことが可能である。

3.5 まとめ

2 進数の有限小数から 10 進数への変換は、常に有限小数になる。しかし例えば、 2^{1023} を 10 進数で表すと、 1 の位の数は 0 ではないので、有効数字は 307 < 1023 log_{10} 2 < 308 より 307 桁になる。従って、実用上は浮動小数点数の表示においても近似が必要である。

ソースが複雑であったため、 printf の中で使われている近似がどのようになされているのか、正確には分からなかった。また、近似無しで浮動小数点数を 10 進数で表示する関数も標準で用意されていない。よって、浮動小数点数を 10 進で正確に、或いは精度保証つきで表示する場合は、自前で関数を用意する必要がある。

4 解答 1-3

これは簡単であったようだ。以下は標準的な解答 (河邊昌彦君) である。問題の式は

\frac{1}{1 + 2x} − \frac{1 − x}{1 + x} = \frac{1 + x − (1 − x)(1 + 2x)}{(1 + 2x)(1 + x)}
= \frac{1 + x − (1 + x − 2x^{2})}{(1 + 2x)(1 + x)}
= \frac{2x^{2}}{(1 + 2x)(1 + x)}

と変形できる。このように変形すれば、同じような数同士の減算が無くなるので、桁落ちを少なくすることができる。

5 解答 1-4

これも簡単であったようだ。以下は標準的な解答 (河邊昌彦君) である。以下は、 Visual C++ 6.0 のマニュアルからの引用である。

volatile として宣言したオブジェクトは、値が常に可変であるため、最適化処理では使われません。システムは volatile オブジェクトの値が要求されると、前回の命令で同一オブジェクトの値を要求されていた場合でも、その時点の値を必ず読み取ります。また、このオブジェクトの値は代入直後に書き込まれます。

即ち、 volatile として変数を宣言すると、プログラムに沿った文の実行以外において参照、変更される可能性がある、ということをコンパイラに示すことになる。プログラムに沿った文の実行以外による参照、変更というのは、他のプロセスやスレッド、割り込みルーチンからのアクセス等がある。このような変数に対しては、プログラムに書いてあるとおりの順番で文を実行しなければ、予想外の結果を引き起こすことになる。従って、コンパイラは、予め値を計算する、実行順を入れ替える、結果をキャッシュする等の最適化を行わなくなる。

ここで、教科書 p.9 のプログラム

     double x1=1.;
     double x10=10.;
     fpsetround(FP_RP);
     x_u=x1/x10;
     fpsetround(FP_RM);
     x_l=x1/x10;

において、 1–2 行目を

     volatile double x1=1.;
     volatile double x10=10.;

と変更したとする。すると、 4 行目で x1/x10 を評価するまでに、 x1, x10 の値が変更される可能性があることになる。このプログラムでは、これらの変数の値が実際に変更されることは無い。しかし、コンパイラは最適化をせず、 fpsetround(FP_RP); の実行より後に x1/x10 を評価するコードを生成する。その結果、丸めモードの指定の後に浮動小数展演算が行われるので、正しい結果が得られるのである。

6 解答 2

これもそれ程難しくなかったようだ。幾つかの例を示そう。

6.1 河邊昌彦君の解答

教科書に沿った部分のみを書く。

2 進数の (規格化) 浮動小数点 a は

a = ± (1 + \frac{d_{2}}{2} + \frac{d_{3}}{2^{2}} + ⋯ + \frac{d_{N}}{2^{N−1}} ) × 2^{e}  (1.1)
d_{i} ∈ {0, 1} (2 ≤ i ≤ N)
e_{min} ≤ e ≤ e_{max}

と表される。また、 (1.1) において

m = ± (1 + \frac{d_{2}}{2} + \frac{d_{3}}{2^{2}} + ⋯ + \frac{d_{N}}{2^{N−1}} )  (1.2)

を符号付仮数と呼び、 e を指数と呼ぶ。このような浮動小数点数の絶対値の最大値は

x_{max} = (1 + \frac{1}{2} + \frac{1}{2^{2}} + ⋯ + \frac{1}{2^{N−1}} ) × 2^{e_{max}}  (1.3)

であり、最小値は

x_{min} = 1 × 2^{e_{min}}  (1.4)

である。倍精度においては、 (1.3), (1.4) は

x_{max} = (2 − 2^{−52}) × 2^{1023} = 1.7976931348623157 ⋯ × 10^{308}  (1.5)
x_{min} = 2^{−1022} = 2.225073858507201 ⋯ × 10^{−308}  (1.6)

となる。零は

+ (1 + \frac{0}{2} + \frac{0}{2^{2}} + ⋯ + \frac{0}{2^{N−1}} ) × 2^{e_{min}}  (1.7)

と表される。仮数部の精度が 53bit であり、

2^{−52} = 2.220446049250313 ⋯ × 10^{−16}  (1.8)

であるから、 e = 0 とすると、

{x | x ∈ R, 1 ≤ x < 2}

の範囲を約 2 × 10^{−16} 刻みで表現できるので、 10 進ではおよそ 15.7 桁の精度がある。

以上は規格化されている浮動小数点数の値の範囲であった。これに対し、非規格化浮動小数点数の絶対値の最小値は、

\frac{1}{2^{53}} × 2^{−1022} = 2^{−1075} = 2.47032822920623 ⋯ × 10^{−324}  (1.9)

である。

6.2 斎藤卓也君の解答 (非常に詳しく精密である)

IEEE Standard for Binary Floating-Point Arithmetic

6.2.1 概要目的この規格の目的は,100%ソフト,100%ハード, 又はソフトとハードの組み合わせによって実現される, 規格を満足する浮動小数点システムの実装にある. この規格によりプログラマやシステムのユーザは, 利用環境が規格を満たしているかどうかという点を判断することができる. 規定事項この規格が規定しているのは

1. 浮動小数点の基本フォーマット, 拡張フォーマット

2. 四則演算・平方根・剰余・比較命令

3. 整数・浮動小数点間のフォーマット変換

4. 異なる浮動小数点フォーマット間の型変換

5. 浮動小数点と 10 進文字列の変換

6. 非数 (NaN) を含めた浮動小数点例外である. 規定外事項以下のことについて, この規格は規定をしていない.

1. 10 進文字列と整数のフォーマット

2. NaN の符号と有効数字 (signi cand) に関する解釈

3. 拡張フォーマットについての 2 進 10 進変換

6.2.2 定義この規格が定義している言葉は以下の通りである.

• biased exponent [げたばき指数] 指数フィールドが非負になるように選択された定数 (bias と呼ぶ) と, 指数の値との和.

• binary oating-point number [2 進浮動小数点数] 符号, 符号付き指数, 有効数字の 3 つの部分からなるビット列. もしそのビット列が数値として存在するならば, その値は 2 の指数乗と有効数字の, 符号つきの積で表す. この規格の中においては, ビット列と, そのビット列が表現している数というものを, 常に区別しているわけではない.

• denormalized number [デノーマル (不正規化数)] 指数が予約された値 (通常は最も小さな値にする) の非ゼロの浮動小数点数のことで, 明示的でも暗示的でも最上位の有効数字のビットが 0 になる. 拡張フォーマットは, 有効数字のけち表現を採用しておらず, 明示的に MSB を 0 にすることで, 正規化されていない値を表現する.

• destination [デスティネーション] 命令の結果のために用意されたロケーションのこと. デスティネーションは, ユーザが明示的に指定しても, システムによって暗示的に用意されていても良い.

• exponent [指数] 2 進浮動小数点の構成要素の 1 つで, 表現している数値を決定する際に,2 の肩に乗る整数. 通常この指数は, 符号付き指数またはバイアスされていない (unbiased) 指数と呼ばれる.

• fraction 有効数字フィールド 2 進表現の有効数字の内, 小数部分を示すもの.

• mode [モード] ユーザがセット, 読み込み, 書き込みが行える値で, 次に続く数値演算命令に実行制御を渡すことのできるもの. デフォルトのモードとは, 明示的に反対を意図したステートメントを, プログラムまたはその仕様書に盛り込むことをしない限り, 反映されると (ユーザが) 仮定してもよいモードのことである. 丸め誤差の方向を制御する丸めモードは, 必ず実装しなければならない (shall). そして実装によっては, 結果の精度を落とすための丸めの精度が要求される. また, オプションとして, 例外を扱うためのトラップイネーブル/ディゼーブルのモードを実装してもよい.

• NaN [Not Number :非数] 浮動小数点フォーマットに基づいてエンコードされるシンボリックな実体である.NaN には,Signaling NaN と Quiet NaN という 2 種類の NaN が存在する.Signaling NaN(s)(複数形で表記, 即ち NaN の表現は一意に決められていない) は, それがオペランドとして現れた時に, 無効演算例外のシグナルを発生する.Quiet NaN(s) は, シグナルを発生しないほとんど全ての数値演算命令で伝播していく.

• result [結果] デスティネーションによってもたらされる結果 (ビット列) で, そのビット列は通常, 数値を表す.

• signi cand [有効数字] 2 進の浮動小数点の構成要素の 1 つで, 暗示的または明示的な小数点よりも左側にある先頭のビット (leading bit) と, 小数点より右側にあるビット列 (fraction) からなる.

• shall この規格の準拠に際して必ず実装しなければならないものについて, そのとき shall という言葉を使用している.

• should 義務ではないが, 実装を強く推奨している時に, この単語が用いられている.

• status ag [ステータスフラグ] セットとクリアができる,2 つの状態をとる値のこと. ユーザはフラグをクリアしたり, コピーしたり, 前の状態に復帰させることができる. それをセットしたとき, ステータスフラグはシステムに依存した情報を含み, あるユーザに対してアクセスができないようにすることができる. 規格の命令は, 以下のフラグのセットを行うことがある.

– 不正確

– アンダーフロー

– オーバーフロー

– 0 除算

– 無効演算

• user [ユーザ] しかしながらこの規格では, 規格で定められたアクセスを行ったり, プロフラミング環境の命令を制御する人やハードウェアやプログラムを特定していない.

6.2.3 フォーマット

規格では,2 つのグループに分けることのできる全 4 種類の浮動小数点フォーマットを定義している. グループには基本フォーマットと拡張フォーマットがあり, 各々について単精度と倍精度を持っている. 規格の実装レベルは, フォーマットのサポートの組み合わせによって判断することができる.

値の集合この節ではフォーマットの範囲内で表現される数値についてのみ規定し, どのようにビット列のエンコードを行うかについては規定しない. そのような数値について, 具体的にどのようなフォーマットを選択するかに関しては,

• 有効数字の桁数である仮数のビット数 p

• 指数の最大値 Emax

• 指数の最小値 Emin

の 3 つの整数で規定することができる. これらの値を含めたパラメータの一覧表を表 1 に示す. そして各々の実体は, 以下のような形で提供しなければならない (shall).

(−1)^{s}2^{E}(b_{0}.b_{1}b_{2} ⋯ b_{p−1})

但し

• s は 0 または 1

• E は E min から E max の間の任意の整数

表 1: 規定されている各フォーマットのパラメータ

パラメータ基本単精度拡張単精度倍精度拡張倍精度
有効数字の桁数 p  24  32  以上 53  64 以上
指数の最大値 Emax  +127  +1023 以上 +1023 +16383 以上
指数の最小値 Emin  − 126  − 1022 以下− 1022 − 16382 以下
指数のバイアス  +127  規定しない +1023 規定しない
指数のビット長  8  11 以上  11  15 以上
全体のビット長  32  43 以上  64  79 以上

• b i は 0 または 1

である. また, 各フォーマットは

• +∞ と −∞ の 2 つの無限大

• 最低 1 種類の Signaling NaN

• 最低 1 種類の Quiet NaN

を持たなければならない.

以上のようなフォーマットの規定では,1 つの数値を表現するのには冗長である. その例として,2^{0}(1.0) = 2^{1}(0.1) = 2^{2}(0.01) が挙げられる. しかし, このような冗長性は拡張フォーマットにおいてのみ起こり得る. そして± 2^{Emin}(0.b_{1}b_{2} ⋯ b_{p − 1}) で表わされる非ゼロの値を不正規化数と呼ぶ. 特定の指数の値は非数 (NaN), 無限大 (∞), ゼロと不正規化数のために予約されている. なお, ゼロには全てのフォーマットについて符号の情報が付加されており,0 除算のようなケースにおいてその効力を発揮する.

基本フォーマット基本単精度と基本倍精度で表現される数は,3 つのフィールドから構成されている.

1. 1 ビットの符号 s

2. げたばき指数 e =E +bias (E はバイアスされていない指数)

3. 有効数字フィールド (fraction)f = .b_{1}b_{2} ⋯ b_{p − 1}

3 つのフィールドの値から, エンコードされる値を表 2, 表 3 に示す.

フォーマットの組み合わせ必ずサポートしなければならないのは, 単精度だけであり (shall), 実装する基本精度 (単精度・倍精度) に対応して, 拡張フォーマット (拡張単精度・拡張倍精度) を実装することを推奨している (should). 何故このような推奨をしているのかについては, 拡張精度を導入する動機の一つとして,10 桁を表示しながらも, 内部では 13 桁の計算を行うような電卓を喩えに挙げて説明することができる. 基本精度に相当する 10 桁に対して,13 桁という拡張精度を導入することで, 電卓はユーザに対して正確で単純な計算モデルを提供することができる. そのため, 単体で拡張精度を導入するのではなく, 基本精度に対応した拡張精度の導入を推奨しているのである.

表 2: 単精度の記憶形式\begin{matrix}s, e, f & 値 \\ s = any, e = 255, f ≠ 0 & NaN \\ s = any, e = 255, f = 0 & (− 1)^{s}∞ \\ s = any, 0 < e < 255, f = any & (− 1)^{s}2^{e − 127}(1.f) \\ s = any, e = 0, f ≠ 0 & (− 1)^{s}2^{e − 126}(0.f) \\ s = any, e = 0, f = 0 & (− 1)^{s}0 \end{matrix}

表 3: 倍精度の記憶形式\begin{matrix}s, e, f & 値 \\ s = any, e = 2047, f ≠ 0 & NaN \\ s = any, e = 2047, f = 0 & (− 1)^{s}∞ \\ s = any, 0 < e < 2047, f = any & (− 1)^{s}2^{e − 1023}(1.f) \\ s = any, e = 0, f ≠ 0 & (− 1)^{s}2^{e − 1022}(0.f) \\ s = any, e = 0, f = 0 & (− 1)^{s}0 \end{matrix}

6.2.4 丸め

丸めとは, まず無限の精度で計算を行ってみて, そして必要ならば, デスティネーションのフォーマットに合わせた形に修正し, 同時に不正確の例外を発生させるものである. 但し,2 進 10 進変換において, 不正確の例外は発生しない. 実際には, この規格で述べられているように無限の精度で計算を行うことは不可能であり, 丸めた結果が常に無限の精度で計算したものと同じになることが保証されるように数ビットの保護ビットを用いてこれに対応する.

丸めモードは, 比較と剰余演算を除いた全ての演算に影響を与え, 結果が 0 になるような加 (減) 算の符号だけでなく, オーバーフロー, アンダーフローのしきい値についても影響する.

最近値への丸めこの規格では, 最近値への丸め (RN :Round to Nearest) をデフォルトの丸めモードとして提供しなければならない (shall). このモードでは表現できる値を, 無限の精度の中間結果に最も近い値にしなければならない (shall). もし, 中間結果に最も近い表現できる値が 2 つあるときは,LSB が 0 になる方を選択しなければならない (shall). しかし, 無限の精度で計算した中間結果が,2^{Emax}(2 − 2^{− p}) を越える場合, 丸めの精度を変更していなければ, 丸めた結果は ∞ となる.

有向丸め最近値への丸めモードだけでなく, 有向丸めモードもユーザに選択権を与えなければならない (shall). 有向丸めモードとは,+∞ 方向への丸め (RP :Round to +∞), −∞ 方向への丸め (RM :Round to −∞),0 方向への丸め (RZ :Round to Zero) の 3 つのモードを指す.

+∞ 方向への丸めは, 無限の精度の中間結果に対して,+∞ の方向に向かってそのフォーマットで最も近い値をとらなければならず, また, 中間結果よりも小さな値になることはない (shall). −∞ 方向への丸めは, 無限の精度の中間結果に対して, −∞ の方向に向かってそのフォーマットで最も近い値をとらなければならず, また, 中間結果よりも大きな値になることはない (shall).0 方向への丸めは, 無限の精度の中間結果に対して,0 の方向に向かってそのフォーマットで最も近い値をとらなければならず, また, 絶対値で比較して, 中間結果よりも大きな値になることはない (shall).

有向丸めにおいて,RP で丸めを行うと, 誤差の上限がわかり,RM は誤差の下限がわかる. また,RZ の丸めを行うと, 誤差の絶対値の上限がわかる.RP と RM の利用価値は, 上限と下限の両方を保証して結果を返す演算 (区間演算) にある. 区間演算を用いると, 誤差の範囲を正確に保証することができる.

丸めの精度通常は, デスティネーションの精度と同じ精度で結果は丸められるものであるが, システムの中には, 倍精度もしくは拡張精度しかデスティネーションが選択できないものがある. そのようなシステムにおいてユーザ (大抵は高水準言語のコンパイラ) は, 単精度で丸められるような仕様にしておかなければならない.

6.2.5 命令

この規格では, 以下の命令を備えなければならない (shall).

• 加減乗除

• 平方根

• 剰余

• 浮動小数点の整数への丸め

• 異なる浮動小数点フォーマットへの変換

• 浮動小数点フォーマットと整数フォーマット間の変換

• 2 進 10 進変換

• 比較

これらの命令は, まず無限の精度で計算を行い, その後デスティネーションのフォーマットに合わせて修正を行わなければならない (shall).

算術演算算術演算に関しては, 全てのフォーマットの 2 オペランドに対して, 加減乗除と剰余演算ができるように実装しなければならない (shall). また,2 つのオペランドのフォーマットが互いに異なるものでも演算ができることが好ましい (should). そして丸めの精度に関わらず, デスティネーションのフォーマットは 2 つのオペランドのフォーマットよりも精度の大きなものにしなければならない. 即ち, デスティネーションの精度は, 入力オペランドの持つ精度よりも粗い精度にしてはいけない. この精度に関して, 入力オペランドが倍精度の 2 次元ノルムである。\sqrt{x^{2} + y^{2}} の計算を例にとって規格外の演算というものを考える. この計算は自乗の和の計算と, 平方根の 2 段階で行われるので,(x^{2} + y^{2}) を拡張倍精度で計算して, それから, 平方根計算について精度を落として倍精度の結果を出したとする. このような演算を内部で行うと, 丸めの問題からシステムによって異なる答えを返す可能性がある. それを防ぐために精度を落とす演算環境は禁止されているのである.

剰余演算 r = (xREMy) において,y =0 のとき,r は丸めモードに関わらず, 数学的な r = x − y × n(n は整数で,x/y の正確な最近値) という関係で定義される. ところで |n − x/y| が正確に\frac{1}{2} のときは, 最近値 n は 2 通りの整数を取り得ることになるが, 規格ではこのような時は,n は偶数とすると定められている. 即ち, 剰余演算は丸めモードや丸めの精度等によらず, 常に正確な一意の値を返す.

平方根演算平方根演算は, 実装されている全てのフォーマットで計算できなければならない (shall). オペランドが正の値の時, 結果の符号は正であるが, − 0 の平方根は, − 0 としなければならない (shall).

浮動小数点フォーマット変換サポートしている全てのフォーマット間について, フォーマットの変換ができなければならない (shall). 低い精度へのフォーマット変換については,4 節で規定した通りに丸めを行わなければならない (shall). また高い精度へのフォーマット変換については, 結果は常に正確である.

浮動小数点 ↔ 整数フォーマット変換サポートされている全ての浮動小数点と整数フォーマットについて変換ができなければならない (shall). 浮動小数点から整数への変換は,4 節で規定された通りの丸めを行わなければならない (shall).7 節で規定している例外が起きない限り, 浮動小数点整数フォーマット変換は正確でなければならない (shall).

浮動小数点の整数への丸め規格では全ての浮動小数点について, 値を整数に丸めて, オペランドと同じ浮動小数点フォーマットに結果を返せなければならない (shall). 丸めの規則は 4 節で規定した通りとする.

２進 ↔ １０進変換最低 1 つの決められたフォーマットでの 10 進文字列とサポートされている全ての 2 進浮動小数点との間で, 変換命令を実装しなければならない (shall). 変換をサポートしなければならない値の範囲は, 表 4,5 の通りとする. 表中の整数 M と N は, ± M × 10^{± N} という値を示すパラメータとする. 但し,M を表の範囲内で, 連続している 0 のビット列 (trailing zeros) を最下位に加えたり, 取り除いたりすることで,N が最小になるように解釈しなければならない. 例えば, 単精度において, 仮に M =1 ,N =107 で表される数値 1 × 10^{107} は,M を調節することで,M =100000000 ,N =99 即ち 100000000 × 10^{99} と解釈できるので, 変換の範囲内として取り扱わなければならない.

表 4: 10 進→ 2 進変換の範囲

フォーマット MAX M MAX N
単精度  10^{9}− 1  99
倍精度  10^{17}− 1  999

表 4,5 と後述する表 6 において, 整数 M が表の範囲を越えたときはオプションで, 単精度なら 9 桁目, 倍精度なら 17 桁目以降の桁を他の数字 (通常は 0) で置き換えるようにすることができる. 表 6 で規定された範囲内において, 変換は 4 節で規定した通りの正確な丸めを行わなければならない.

もし, 表 6 の範囲外の変換を行う際, 丸めモードが最近値への丸めモードになっているときは変換による誤差が 0.47ulp(Unit in the Last Place :デスティネーションの LSB を 1 単位とした値) を

表 5: 2 進→ 10 進変換の範囲

フォーマット MAX M MAX N
単精度  10^{9}− 1  53
倍精度  10^{17}− 1  340

越えてはならない (shall). 丸めモードが有向丸めモードのときは, 正確な符号と 1.47ulp 以下の誤差を保証しなければならない (shall).

表 6: 10 進 ↔2 進変換における正確な丸めの保証範囲

フォーマット MAX M MAX N
単精度  10^{9}− 1  13
倍精度  10^{17}− 1  27

変換命令は,2 進浮動小数点の値が増加したにも関わらず, それらを 10 進に変換を行った際, 値が減少することは許されない (shall). 同様に, 増加していく 10 進の値を 2 進に変換する際, その 2 進の値が減少することは許されない (shall). 即ち,2 進で保証されている大小関係は,10 進に変換しても保たれなければならない.

最近値への丸めモードで, 値が表 4,5 の範囲内にある 2 進浮動小数点を 10 進に変換して, それをまた再び 2 進数に戻した時, 変換後の 2 進数から元の 2 進数を一意に識別できなければならない (shall). また, 表 4,5 の範囲を言い換えれば, この規則の適用範囲は単精度なら 10 進数で 9 桁, 倍精度なら 10 進数で 17 桁以内の数値ということになる.

もし,10 進から 2 進への変換でオーバーフロー/アンダーフローが起きた場合,7 節で定められた通りの例外処理を行う必要がある.2 進から 10 進の変換の際, オーバーフロー, アンダーフロー, 非数 (NaN), 無限大 (∞) と判断される場合は, それらを示す適切な文字列をユーザに提供することが好ましい (should). また,10 進文字列の NaN の表現は, この規格では規定しない.

一貫性を保つために,2 進 ↔10 進変換のプロシージャは, 実際に変換を行うのが, インタプリタ・コンパイラ・アセンブラ等の言語変換なのかランタイム・対話的入出力等のプログラムの実行なのかによらず, 常に同じ結果を返すことが好ましい (should).

比較比較命令は, サポートされている全ての浮動小数点について行えなければならない (shall). そしてそれは異なるフォーマット間の比較についても適用される. 比較は常に正確で, オーバーフローやアンダーフローを引き起こすことはない.4 つの排他的な関係にある条件 less than,equal,greater than,unordered が設定できる. 最後のケースである unordered は, 少なくとも 1 つのオペランドが NaN のときの比較の際に発生する条件であり,NaN はどのような値との比較も unordered でなければならない (shall). これを利用した例として, ある値 x が NaN かどうか調べるときは,x と x を比較して unordered なら NaN であると判断できる. 比較に際しては, ゼロの符号は無視, 即ち − 0 = +0 でなければならない (shall).

− 0 = +0 でなければならない (shall). 比較の結果は, 実装者の二者択一のオプションとして規定されている. 一つは, 上記の 4 つの条件を特定することができる条件コードを結果とするものであり, もう一つは, 特定した条件 (predicate) に対する比較に true/false を返す方法である. また,true/false 型の場合, 比較結果が unordered の時は, 表 7 の条件 (最終列の YES/NO) に基いて, 後述する無効演算の例外を出さなければならない. 表中の特定条件の項目は,ad hoc 表現に加えて FORTRAN ライクにも書かれており, 関係の項目の T/F は,True/False を意味し, 比較不能の条件は?で表記しており,4 つの T/F で (FORTRAN ライク表現と等価な)1 つの状態を表していると考えることができる.

比較命令としてサポートしなければならないのは, 表 7 の内, 最初の 6 つであり (shall), 更に 7 つ目の条件 (unordered ) の実装を推奨している (should). また, これに対応する NOT() の命令についても同様である.

表 7: 比較命令の対応

特定条件 (in FORTRAN) 関係 (GT,LT,EQ,UN) 例外発生?
=(.EQ.)  (F,F,T,F)  NO
? <> (.NE.)  (T,T,F,T)  NO
> (.GT.)  (T,F,F,F)  YES
>= (.GE.)  (T,F,T,F)  YES
(.LT.)  (F,T,F,F)  YES
=(.LE.)  (F,T,T,F)  YES
?(unordered)  (F,F,F,T)  NO
<> (.LG.)  (T,T,F,F)  YES
=¿(.LEG.)  (T,T,T,F)  YES
? > (.UG.)  (T,F,F,T)  NO
? >= (.UGE.)  (T,F,T,T)  NO
?(.UL.)  (F,T,F,T)  NO
? = (.ULE.)  (F,T,T,T)  NO
? = (.UE.)  (F,F,T,T)  NO
NOT (>)  (F,T,T,T)  YES
NOT (>=)  (F,T,F,T)  YES
NOT ()  (T,F,T,T)  YES
NOT (=)  (T,F,F,T)  YES
NOT (?)  (T,T,T,F)  NO
NOT (<>)  (F,F,T,T)  YES
NOT (=>)  (F,F,F,T)  YES
NOT (? >)  (F,T,T,F)  NO
NOT (? >=)  (F,T,F,F)  NO
NOT (?)  (T,F,T,F)  NO
NOT (? =)  (T,F,F,F)  NO
NOT (? =)  (T,T,F,F)  NO

6.2.6 無限大, 非数, 符号付きゼロ

∞に関する演算無限大 (∞,INF) は, −∞(あらゆる有限数)+∞ として解釈しなければならない (shall). ∞ に対する演算は, 常に正確であり, かつ 7 節で規定されている演算を除いて, 例外のシグナルを発生させることはない (shall). ∞ に付随する例外は以下の通りである.

1. ∞がオーバーフロー,0 除算によって生成され, かつトラップがディゼーブルであるとき

2. ∞が無効演算である場合

NaN に対する命令 NaN については, 全ての命令に対して,Signaling NaN と Quiet NaN の異なる 2 種類の NaN をサポートしなければならない (shall).Signaling NaN は, 初期化されていない変数や, 規格で規定されていない, 虚数や巨大な範囲を扱うといった規格の拡張の際に用いられる.Quiet NaN は, 無効若しくは利用不可能なデータや演算結果を, さかのぼって判定する診断情報として利用することができる (should). それ以外の Quiet NaN の利用に関しては, 実装者の裁量に任せられる. 診断情報の伝播において要求されることは,NaN(s) に含まれている情報を, 数値演算命令と浮動小数点フォーマット変換命令の際に保存しておくことである (should).

Signaling NaN(s) は,5 節で挙げた全ての命令について, 後述する無効演算例外のシグナルを送るオペランドとして予約されている (shall). フォーマットの変更なく,Signaling NaN をコピーしたときに, 無効演算例外を発生させるかどうかは, 実装者に選択権がある.

Signaling NaN を含んだ全ての演算と, 無効演算について, もしトラップが起こらないで, かつ浮動小数点の結果が返されるべきものであるならば, 結果として Quiet NaN を返さなければならない (shall).

入力に 1 つないし 2 つの NaN を含む全ての演算は, そのどちらも Signaling NaN でないときは, 例外のシグナルは発してはならない. そして, 浮動小数点の結果が返されるべきものであるならば, 結果として Quiet NaN を返さなければならない (shall). そして結果として返すその Quiet NaN は, 入力のうちのどちらかの NaN を, そのまま返すことが望ましい (should). しかしフォーマット変換の場合は, 全く同じ NaN を伝播させることはできないことに注意しなければならない.Quiet NaNs は, 浮動小数点フォーマットの結果を返さない命令における Signaling NaNs のふるまいと, ほぼ同じ効果を持っている (その命令とは即ち, 既に述べられている,NaN を伴わない比較, フォーマットの変換のことである).

符号ビット規格では,NaN の符号を解釈することはない. その NaN 以外のケースにおいて

• 積と商の符号は,2 つのオペランドの符号の xor をとった値である

• 減算 x − y を x + (− y) という加算とみなすと, 被加算数の符号が 2 つとも同じ符号であるならば, 和の符号もそれと同じ符号をとらなければならない

• 浮動小数点を整数に丸める命令の結果の符号は, オペランドの符号と等しい

• 異符号同士の和 (または同符号同士の差) が丁度 0 になったとき, その和 (差) は

– 丸めモードが RM( −∞ 方向への丸め) のときは− 0

– それ以外のモードの場合は+0 を返す • しかし,x = ±0 のときに x + x または x − (− x) を行ったら, その x の符号が保持されたま

まそれが結果となる, 即ち, 結果として x をそのまま返す上の 2 点をまとめると, 表 8 のようになる

• \sqrt{− 0} = − 0 であることを除いて, 全ての平方根演算は, 正の符号を返す

表 8: 結果が 0 となる加減算演算

計算丸めモード結果

(+0)ADD(+0)

(+0)SUB(− 0) (don ’ t care) +0

( − 0)ADD( − 0)

( − 0)SUB(+0) (don ’ t care) -0

(+0)ADD( − 0) (+0)SUB(+0) ( − 0)ADD(+0) ( − 0)SUB( − 0) x SUB x

x ADD ( − x ) RN,RZ,RP +0

(+0)ADD( − 0) (+0)SUB(+0) ( − 0)ADD(+0) ( − 0)SUB( − 0) x SUB x

x ADD ( − x )  RM  -0

6.2.7 例外

規定された 5 種類の例外のいずれかを検出したときに, シグナルを発生しなければならない (shall). シグナルとは, ステータスフラグを立てる, または, トラップを発生すること, 若しくはその両方のことである. 各々の例外のトラップは,8 節で規定しているように, ユーザの制御下におくことを推奨する (should). 例外に対するデフォルトの応答は, トラップ無しに処理を続行させることである (shall). 規格では, 結果について, トラップと非トラップの両方を規定している. そしてトラップがイネーブルかどうかで結果が変化するものがある (この代表的なものが, アンダーフローに対するふるまいである).

各々の例外について, 対応したトラップが発生しないとき, 代わりに値がセットされるようなステータスフラグを提供しなければならない (shall). このステータスフラグは, ユーザの要求によってのみリセットされる (スティッキーである). ユーザはステータスフラグを独立にテスト, 入れ替えができ, なおかつ 5 つのフラグを一度にセーブ, リストアできなければならない (shall).

複数同時に起こり得る例外は, 不正確とオーバーフローの例外と, 不正確とアンダーフローの例外の 2 つだけである (つまりそれ以外の例外は 1 つずつしか起こらない).

無効演算規格では, 命令を遂行することが意味が無いようなオペランドを持つときに, 無効演算例外が発生する. その結果は, トラップがなく例外が発生したとき,Quiet NaN を浮動小数点に対応したフォーマットでデスティネーションに提供しなければならない (shall). 無効演算の種類は以下の通りである.

• Signaling NaN に関する全ての演算

• (+∞) − (+∞) のような無効な ±∞ の加減算

• 0 × ∞

• 0 /0 と ∞/∞

• 剰余演算 xREM y において,x が ∞ または y が 0 であるとき

• 0 より小さな値の平方根演算

• 2 進浮動小数点から, 整数または 10 進のフォーマットに変換する際, オーバーフロー・ ∞ ・ NaN によって, 正しい変換が行うことができないとき

• 比較不能 (unordered ) な 2 数について,?(unordered ) を含まない条件 (¡¿等) で比較命令を行ったとき

0 除算除数が 0 で, 被除数が有限の非ゼロの数であるならば,0 除算例外をシグナルを発しなければならない (shall). そして, トラップが起こらないならば, 結果として適切な符号付きの∞を返さなければならない (shall).

オーバーフローオーバーフローの例外シグナルは, 指数を限りなくとれると仮定して丸められた浮動小数点の結果が, 絶対値で比較して, デスティネーションのフォーマットの最大の有限の数を越えた場合に発生する (shall). そしてその結果は, トラップが起こらないなら, 丸めモードと, 中間結果の符号によって決定される.

1. 丸めモードが RN なら, オーバーフローのとき, 中間結果の符号付き∞を返す.

2. 丸めモードが RZ なら, オーバーフローのとき, フォーマットで表現できる最大値を, 中間結果の符号をつけて返す.

3. 丸めモードが RM なら, 正のオーバーフローのとき, フォーマットで表現できる最大値を返し, 負のオーバーフローのとき, − ∞ を返す.

4. 丸めモードが RP なら, 負のオーバーフローのとき, フォーマットで表現できる負の最大値を返し, 正のオーバーフローのとき,+∞ を返す.

以上をまとめると, 表 9 のようになる.

フォーマット等の変換命令を除く全ての命令について, トラップを発生させるオーバーフローは, 無限の精度で計算した中間結果を,2α で割り, それを丸めたものをハンドラに渡さなければならない (shall). ここで α とは, 指数フィールドが 8 ビットの単精度なら 192 ,11 ビットの倍精度なら 1536,n ビットの拡張精度なら 3 × 2^{n − 2} である. つまり, 正確に計算した中間結果の指数フィールドから,3 × 2^{n − 2} を引いて結果を返し, かつフォーマットに合うように仮数を正規化して返す.2 進表 9: オーバーフロー時における丸めモードと結果の関係

丸めモード中間結果の符合結果
RN  +  +∞
RN  -  −∞
RZ  +  +(表現しうる最大値)
RZ  -  -(表現しうる最大値)
RP  +  +∞
RP  -  -(表現しうる最大値)
RM  +  +(表現しうる最大値)
RM  -  −∞

浮動小数点からのフォーマット変換において, トラップ型のオーバーフローが起きた場合には, できるだけデスティネーションの精度で丸めを行いながらも, 指数のバイアスを調整して, その 2 進のフォーマットと同じかそれよりも広範囲のフォーマットを結果としてハンドラに送らなければならない (shall).10 進から 2 進への変換でトラップ型オーバーフローが起きたときは, できるだけデスティネーションの精度で丸めを行いながらも, 指数のバイアスを調整して, サポートされている最も大きなフォーマットを結果としてハンドラに送らなければならない (shall). しかし結果が表現できる範囲から極めて離れていて調整することができない場合は,quiet NaN を代わりに送らなければならない (shall).

アンダーフローアンダーフローには 2 つの相関のあるイベントが関係している.1 つは ±2^{Emin} の範囲内にある, 小さな非ゼロの結果の創出で, これは, 余りにも小さすぎて, 例えば割り算でオーバーフローを引き起こすように, 後で他の例外を引き起こす可能性があるというものである. もう 1 つはデノーマルによる, 非常に小さな数の近似のために, 極端に精度が悪化してしまうことである. 実装者はどのようにそれらのイベントを検出するかを, 与えられた方法の中から選択可能であるが, 全ての命令について一貫した方法でそれらのイベントを検出しなければならない (shall). その検出法は 4 種類に分けることができる. 最初の 2 つが値の小ささ (Tininess) についての判定基準で, 後半 2 つが精度落ち (Loss of Accuracy) についての判断基準である.

1. 丸め後検出 (After rounding)

指数の範囲が無限であるかのようにして (有効数字は有限=丸めて) 計算した非ゼロの結果が, 厳密に ±2^{Emin} の間に存在するとき

2. 丸め前検出 (Before rounding)

指数と精度の両方の範囲が無限であるかのようにして (=丸めないで) 計算した非ゼロの結果が, 厳密に ±2^{Emin} の間に存在するとき

3. デノーマルロスによる検出 (A denormalization loss) 計算した結果と, 指数の範囲が無限であると仮定した結果が異なっているとき, 即ち, ある値を表現するのに, 正規化数と非正規化数とで異なった場合に検出される

4. 不正確による検出 (An inexact result)

計算した結果と, 指数と仮数の両方の範囲が無限であると仮定した結果が異なっているとき (この条件を不正確として後述している. つまり, 不正確であるならば検出するという方法である)

丸め前検出と, 丸め後検出で結果が異なる例としては,2^{Emin} よりも僅かに小さいが, 丸めるとちょうど 2^{Emin} になってしまう数がある. このような数は, 丸め前検出のときは Tininess が検出されず, 丸め後検出では Tininess が検出される. アンダーフロートラップが実装されていないとき, 又はイネーブルされていないとき (これがデフォルトである), 値が小さく (Tininess) かつ精度が失なわれたとき (Loss of Accuracy) のみ, アンダーフローはシグナルを発生させる (アンダーフローフラグを立てる). つまり計算結果がデノーマルになっても, その結果が丸められておらず正確ならば, アンダーフローフラグを立ててはならない (shall). Tininess と Loss of Accuracy を検出するための方法は, 結果に影響を与えてはいけない. また, その結果とは, ゼロかデノーマルか ±2^{Emin} のはずである. アンダーフロートラップが実装され, かつイネーブルの時,Tininess が検出されていたなら,Loss of Accuracy の検出の有無に関わらすアンダーフローのシグナルを発生しなければならない (shall). フォーマット等の変換命令を除く全ての命令について, トラップを発生させるアンダーフローは, 無限の精度で計算した中間結果に 2α をかけ, それを丸めたものをハンドラに渡さなければならない (shall). ここで α とは, 指数フィールドが 8 ビットの単精度なら 192 ,11 ビットの倍精度なら 1536,n ビットの拡張精度なら 3 × 2^{n − 2} である. つまり, 正確に計算した中間結果の指数フィールドから,3 × 2^{n − 2} を足して結果を返し, しかもフォーマットに合うように仮数を正規化して返さなければならない. 変換命令において, トラップを起こすアンダーフローは, トラップを起こすオーバーフローの時の動作とアナロジーな動作をしなければならない (shall).

不正確命令の結果を丸めたものが正確ではないとき, 又はオーバーフローのトラップを起こさない状態で結果がオーバーフローしたとき, 不正確例外をシグナルしなければならない (shall). また, 丸められた若しくはオーバーフローした結果をデスティネーションに送られなければならない (shall). もし, 不正確のトラップが発生するならば, それをトラップハンドラに渡さなければならない (shall).

6.2.8 トラップ

ユーザは,5 つのうち任意の例外が起きたとき, 決められたハンドラに処理を渡すという形でトラップを要求できる (should). また, 存在しているハンドラについて, ディゼーブル, 保存, リストアを要求できる (should). ユーザは, 意図した例外のトラップハンドラがイネーブルされているかどうかを判定することができる (should). トラップがディゼーブルで, 例外シグナルが発生したときは,7 節で規定した方法に従わなければならない.(shall). トラップがイネーブルの状態で, 例外のシグナルが発生したときは, そのプログラムの実行はサスペンドされ, ユーザによって前もって用意されているトラップハンドラを実行しなければならない (shall), そして,7 節で規定された例外の規則に従って結果を返さなければならない (shall).

トラップハンドラトラップハンドラは, 例外を発生させた命令に代わって結果の値を返すサブルーチンの能力を持たせることを推奨している (should). そしてトラップハンドラが結果を返さないならば, その値は未定義とする. 同様に, トラップがイネーブルの状態でシグナルが発生したとき, 例外に対応した各々のフラグは, トラップハンドラによってセット又はリセットされていないならば未定義とする.

トラップが発生したとき, ハンドラは以下のことについて判断できる (should).

1. どの例外が起きたか (例外の種類)

2. どの命令で起きたか (命令の種類)

3. 結果のフォーマット

4. オーバーフロー, アンダーフロー, 不正確の例外の場合は, デスティネーションのフォーマットに合わないという情報を含めた, 正確に丸められた結果

5. 無効演算若しくは 0 除算による例外の場合は, オペランドの値

トラップの優先また, トラップがイネーブルであるならば, オーバーフローとアンダーフローは, 同時に起きた不正確のトラップに対して優先権を持つ. これらの例外以外は同時に発生することがないので, これら 2 つのトラップのみ優先権を与えておけばよい.

6.2.9 付録:推奨関数

IEEE754 規格には付録として, 異なるシステム間でのプログラムの移植性の観点から実装を推奨する 10 個の関数が挙げられている. しかし同時に, この付録は規格の一部ではないことも主張している. 以下に示した関数の中で,(1),(2),(6) のようにフォーマットの変換が行われない関数は,signaling NaN に対して無効演算の例外を発生するかどうかについては実装者の権限に任せられる.

1. Copysign(x,y )

y の符号を,x にコピーして返す. 例えば,Copysign(a,1 .0) を実行すると,a の絶対値 abs(a ) が返ってくる.

2. Reversesign(x)

入力 x の符号を反転させて返す. この関数は 0 − x とは異なる. 具体的には,x = ±0, x = NaN のときに違いが見られる.

• x = +0 のとき,0 − x は, 丸めモードが RM 以外のとき,+0 を返す

• x = NaN のとき,0 − x は,NaN を返す

3. Scalb(y,N )

y × 2^{N} を返す関数.N は整数で, 具体的に 2^{N} という値を計算することなく, 直接浮動小数点を生成する.

4. Logb(x)

x のげたばきでない指数 (符号付き整数) を返す. 但し

• Logb(NaN)=NaN

• Logb( ∞)=+ ∞

• Logb(0)= −∞ (更に 0 除算例外を発生する)

である. x が正の有限数のとき,Scalb [x, − Logb(x )] は厳密に 0 と 2 の間にあり, 特に 1 より小さくなるのは不正規化数のときだけである

5. Nextafter(x,y )

表現可能な数値を並べたとき,y への方向に対して,x の隣にある数値を返す.

• x =y ならば, いかなる例外も起こすことなく x を返す.

• x,y どちらかが Quiet NaN なら, 入力と全く同じ NaN かそれ以外の NaNs を返す

• x が有限でも,Nextafter(x,y ) が ∞ のときは, オーバーフローと不正確の例外を出す

• Nextafter(x,y ) が, 厳密に ±2^{Emin} の中にあるなら, アンダーフローと不正確の例外を出す

6. Finite(x )

−∞ < x < +∞なら TRUE を, そうでなければ FALSE を返す.

7. Isnan(x )

x が NaN かどうかを調べる. 実際には,x と x を比較してその結果,equal にならなければ NaN であると判断できる.

8. GraterOrLess(x,y ) x > y または x < y のときのみ TRUE を返す. 比較不能 (unordered ) という条件があるので, 「x と y は同じではない (x = y)」と区別しなければいけない

9. Unordered(x,y ) 比較不能のときに TRUE を返す.

10. Class(x ) x を 10 種類のクラスに分類する.

• Signaling NaN

• Quiet NaN

• −∞

• 非ゼロ, 負の正規化数

• 負のデノーマル

• -0

• +0

• 正のデノーマル

• 非ゼロ, 正の正規化数

• +∞

この関数は, 例え入力が Signaling NaN であったとしても, 決して例外は発生しない.

7 解答 3

これは難しかったようだ。私のパソコンは Pentium III 750MHz であるが，これには規格違反がある。 4,5 端末室のパソコンも PentiumIII 800MHz で規格違反があるようだ。第一端末室と第６端末室は規格違反がないとでるようであった。

7.1 河邊昌彦君の解答

7.1.1 概要

最近値への丸めを規定している、 [2] の 4.1 節を引用する。

最近値への丸め

この規格では, 最近値への丸め (RN : Round to Nearest) をデフォルトの丸めモードとして提供しなければならない (shall). このモードでは表現できる値を, 無限の精度の中間結果に最も近い値にしなければならない (shall). もし, 中間結果に最も近い表現できる値が 2 つあるときは,LSB が 0 になる方を選択しなければならない (shall). しかし, 無限の精度で計算した中間結果が,2^{E_{max}}(2 − 2^{−p}) を越える場合, 丸めの精度を変更していなければ, 丸めた結果は ∞ となる.

この規定の内、 “このモードでは表現できる値を, 無限の精度の中間結果に最も近い値にしなければならない (shall)” という部分について、検証する。

具体的には、まず、隣接する二つの浮動小数点数をとり、 d_{1}, d_{2}(d_{1} < d_{2}) とする。次に、 (d_{2}−d_{1})/2 より僅かに大きな浮動小数点数をとり、 d_{3} とする。そして、 d_{1} + d_{3} を計算するが、この中間結果は浮動小数点数の精度を越えているため、丸めが行われる。その結果は、正しく最近値への丸めが行われているならば、 d_{2} になるはずである。実際にそのようになっているかを、精度を変えながら試してみる。

7.1.2 プログラムの解説

このテストを行うのが、 7.1.4 のプログラムである。例えば、倍精度浮動小数点数の場合は、

fp_{1} = 1.5 = 1.1_{2} × 2^{0}

ととると、これに隣接する浮動小数点数は

fp_{2} = 1. \underbrace{1000000000000000000000000000000000000000000000000001}_{52bits}_{2} × 2^{0}

である。その差は

fp_{3} = fp_{2} − fp_{1} = 1.0 × 2^{−52}

であり、その半分は

fp_{4} = \frac{fp_{3}}{2} = 1.0 × 2^{−53}

である。ここで、例えば、

fp_{5} = 1.0 × 2^{−52−53} = 1.0 × 2^{−105}

とし、

fp_{6} = fp_{4} + fp_{5}
= 1. \underbrace{0000000000000000000000000000000000000000000000000001}_{52bits}_{2} × 2^{−53}

とする。そして、 fp_{1} と fp_{6} の和を計算する。

fp_{1} + fp_{6} = 1.1_{2} + 1.0000000000000000000000000000000000000000000000000001_{2} × 2^{−53}
= 1.1000 ⋯ 0000_{2}
+ 0.\underbrace{0000 ⋯ 0000}_{52bits}1\underbrace{0000 ⋯ 0001}_{52bits}_{2}
= 1.\underbrace{1000 ⋯ 0000}_{52bits}1\underbrace{0000 ⋯ 0001}_{52bits}_{2}

従って、最近値への丸めを行った結果は、

fp_{7} = fp_{1}□+ fp_{6}
= □ (fp_{1} + fp_{6})
= □ (1.\underbrace{1000 ⋯ 0000}_{52bits}1\underbrace{0000 ⋯ 0001}_{52bits}_{2})
= 1.\underbrace{1000 ⋯ 0001}_{52bits}_{2}

となるので、

fp_{7} = fp_{2}

となる。

7.1.3 実行結果

PentiumIII 733MHz 上でプログラムを実行した結果を 7.1.4 に載せた。それによると、倍精度浮動小数点数を 64bit の計算精度で計算すると、正しく丸めが行われないことがあるようである。詳しく見ると、

   1.1000000000000000000000000000000000000000000000000000e+0
 + 1.0000000001000000000000000000000000000000000000000000e-53
 = 1.1000000000000000000000000000000000000000000000000001e+0

では成功しているが、

   1.1000000000000000000000000000000000000000000000000000e+0
 + 1.0000000000100000000000000000000000000000000000000000e-53
 = 1.1000000000000000000000000000000000000000000000000000e+0

では失敗している。加算する数の最も下位にある 1 のビットは、上の例では 2^{−63} であり、下の例では 2^{−64} である。また、 Pentium プロセッサの拡張倍精度は、仮数部は 64bit の精度である。即ち、倍精度を 64bit の精度で計算すると、拡張倍精度の精度でしか最近値への丸めが行われないようである。しかし、倍精度を 53bit の精度で計算した場合や、拡張倍精度を 64bit の精度で計算した場合は、このようなことは起こらない。よって、内部 64bit→ デスティネーション 53bit の変換が行われることが原因と考えられる。

以上をまとめると、少なくとも PentiumIII プロセッサにおいては、 nearest mode の規格に従っていない部分がある、と結論付けられる。

7.1.4 プログラム

Test nearest mode

ソースコード問題 3 の nearest mode をチェックするプログラムである。 Visual C++又は Borland C++用であるが、 80bit の拡張倍精度をチェックする場合は Borland C++でコンパイルする必要がある (Visual C++の long double は double にマップされているため)。

#include <float.h>
#include <math.h>
#include <typeinfo.h>
#include <iostream>
#include <ostream>
#include <iomanip>
#include <limits>
using std::cout;
using std::endl;
///////////////////////////////////////////////////////////////////////////////
#ifdef __BORLANDC__
void    SetUpRoundMode( void ){ _controlfp( RC_UP, MCW_RC ); }
void    SetDownRoundMode( void ){ _controlfp( RC_DOWN, MCW_RC ); }
void    SetNearRoundMode( void ){ _controlfp( RC_NEAR, MCW_RC ); }
void    SetChopRoundMode( void ){ _controlfp( RC_CHOP, MCW_RC ); }
void    Set64bitPrecision( void ){ _controlfp( PC_64, MCW_PC ); }
void    Set53bitPrecision( void ){ _controlfp( PC_53, MCW_PC ); }
void    Set24bitPrecision( void ){ _controlfp( PC_24, MCW_PC ); }
#endif
#ifdef _MSC_VER
void    SetUpRoundMode( void ){ _controlfp( _RC_UP, _MCW_RC ); }
void    SetDownRoundMode( void ){ _controlfp( _RC_DOWN, _MCW_RC ); }
void    SetNearRoundMode( void ){ _controlfp( _RC_NEAR, _MCW_RC ); }
void    SetChopRoundMode( void ){ _controlfp( _RC_CHOP, _MCW_RC ); }
void    Set64bitPrecision( void ){ _controlfp( _PC_64, _MCW_PC ); }
void    Set53bitPrecision( void ){ _controlfp( _PC_53, _MCW_PC ); }
void    Set24bitPrecision( void ){ _controlfp( _PC_24, _MCW_PC ); }
#endif
///////////////////////////////////////////////////////////////////////////////
typedef unsigned char       BYTE;
typedef unsigned short      WORD;
typedef unsigned int        DWORD;
typedef unsigned __int64    QWORD;
///////////////////////////////////////////////////////////////////////////////
template< class T >

class FloatingPointTraits
{
public:
struct Number{};
};
template<>
class FloatingPointTraits< float >
{
public:
enum { ExponentBias = 127, FractionBits = 23, ExponentBits = 8 };
struct Number
{
    union
    {
        float   fp;
        struct
        {
            DWORD       fraction    : FractionBits;
            DWORD       exponent    : ExponentBits;
            DWORD       sign        : 1;
        };
    };
};
static bool IsNaN( Number fp ){ return ( fp.exponent == 255 && fp.fraction != 0 ); }
static bool IsInfinity( Number fp ){ return ( fp.exponent == 255 && fp.fraction == 0 ); }
static bool IsNormalized( Number fp ){ return ( 0 < fp.exponent && fp.exponent < 255 ); }
static bool IsDenormalized( Number fp ){ return ( fp.exponent == 0 && fp.fraction != 0 ); }
static bool IsZero( Number fp ) { return ( fp.exponent == 0 && fp.fraction == 0 ); }
};
template<>
class FloatingPointTraits< double >
{
public:
enum { ExponentBias = 1023, FractionBits = 52, ExponentBits = 11 };
struct Number
{
    union
    {
        double  fp;
        struct
        {
            QWORD       fraction    : FractionBits;
            QWORD       exponent    : ExponentBits;
            QWORD       sign        : 1;
        };
    };
};
static bool IsNaN( Number fp ){ return ( fp.exponent == 2047 && fp.fraction != 0 ); }
static bool IsInfinity( Number fp ){ return ( fp.exponent == 2047 && fp.fraction == 0 ); }
static bool IsNormalized( Number fp ){ return ( 0 < fp.exponent && fp.exponent < 2047 ); }
static bool IsDenormalized( Number fp ){ return ( fp.exponent == 0 && fp.fraction != 0 ); }
static bool IsZero( Number fp ){ return ( fp.exponent == 0 && fp.fraction == 0 ); }

    static double   GetNextAfter( double fp1, double fp2 ){ return _nextafter( fp1, fp2 ); }
};
#ifdef __BORLANDC__
template<>
class FloatingPointTraits< long double >
{
public:
    enum { ExponentBias = 0x3fff, FractionBits = 63, ExponentBits = 15 };
    struct Number
    {
        union
        {
            long double fp;
            struct
            {
                QWORD       fraction    : FractionBits;
                QWORD       fractionMSB : 1;
                QWORD       exponent    : ExponentBits;
                QWORD       sign        : 1;
            };
        };
    };
    static bool IsNaN( Number fp ){ return ( fp.exponent == 0x7fff && fp.fraction != 0 ); }
    static bool IsInfinity( Number fp ){ return ( fp.exponent == 0x7fff && fp.fraction == 0 ); }
    static bool IsNormalized( Number fp ){ return ( fp.fractionMSB == 1 ); }
    static bool IsDenormalized( Number fp ){ return ( fp.fractionMSB == 0 ); }
    static bool IsZero( Number fp ){ return ( fp.exponent == 0 && fp.fraction == 0 ); }
};
#endif
#ifdef _MSC_VER
template<>
class FloatingPointTraits< long double > : public FloatingPointTraits< double >{};
#endif
///////////////////////////////////////////////////////////////////////////////
template< class T >
class FloatingPoint : public FloatingPointTraits< T >::Number
{
public:
    typedef FloatingPoint< T >          Self;
    typedef FloatingPointTraits< T >    Traits;
    typedef std::numeric_limits< T >    Limits;
    typedef T   ValueType;
    enum
    {
        ExponentBias = Traits::ExponentBias,
        FractionBits = Traits::FractionBits,

        ExponentBits = Traits::ExponentBits
    };
    FloatingPoint( void ){ fp = T(); }
    FloatingPoint( T v ){ fp = v; }
    FloatingPoint( QWORD s, QWORD f, QWORD e ){ sign = s; fraction = f; exponent = e; }
    bool    IsNaN( void ) const { return Traits::IsNaN( *this ); }
    bool    IsInfinity( void ) const { return Traits::IsInfinity( *this ); }
    bool    IsNormalized( void ) const { return Traits::IsNormalized( *this ); }
    bool    IsDenormalized( void ) const { return Traits::IsDenormalized( *this ); }
    bool    IsZero( void ) const { return Traits::IsZero( *this ); }
    int     GetSign( void ) const { return ( sign ? -1 : 1 ); }
    Self    GetNext( void ) const { return Traits::GetNextAfter( fp, Limits::infinity() ); }
    Self    GetBefore( void ) const { return Traits::GetNextAfter( fp, -Limits::infinity() ); }
    int GetUnbiasedExponent( void ) const
    {
        return ( IsNormalized() ? ( int )exponent - ExponentBias : -ExponentBias + 1 );
    }
};
///////////////////////////////////////////////////////////////////////////////
template< class T >
FloatingPoint< T >  operator+( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp + fp2.fp ); }
template< class T >
FloatingPoint< T >  operator-( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp - fp2.fp ); }
template< class T >
FloatingPoint< T >  operator*( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp * fp2.fp ); }
template< class T >
FloatingPoint< T >  operator/( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp / fp2.fp ); }
template< class T >
FloatingPoint< T >  sqrt( FloatingPoint< T > fp )
{ return ( sqrt( fp1.fp ) ); }
template< class T >
bool    operator==( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp == fp2.fp ); }
template< class T >
bool    operator!=( FloatingPoint< T > fp1, FloatingPoint< T > fp2 )
{ return ( fp1.fp != fp2.fp ); }
///////////////////////////////////////////////////////////////////////////////
typedef FloatingPoint< float >          Float;
typedef FloatingPoint< double >         Double;
typedef FloatingPoint< long double >    LongDouble;

///////////////////////////////////////////////////////////////////////////////
template< class T >
std::ostream    &operator<<( std::ostream &ostr, FloatingPoint< T > fp )
{
if ( fp.GetSign() < 0 )
    ostr << "-";
if ( fp.IsZero() )
    ostr << "0";
else if ( fp.IsInfinity() )
    ostr << "Infinity";
else if ( fp.IsNaN() )
    ostr << "NaN";
else
{
    ostr << ( fp.IsDenormalized() ? "0." : "1." );
    for ( int i = fp.FractionBits - 1 ; i >= 0 ; --i )
        ostr << ( int )( ( fp.fraction >> i ) & 1 );
    ostr << "e" << std::setiosflags( std::ios_base::showpos )
        << fp.GetUnbiasedExponent()
        << std::resetiosflags( std::ios_base::showpos );
}
return ostr;
}
///////////////////////////////////////////////////////////////////////////////
template< class T >
T   Test( int bits, int prec )
{
if ( prec == 24 )
    Set24bitPrecision();
else if ( prec == 53 )
    Set53bitPrecision();
else
    Set64bitPrecision(), prec = 64;

T   fp1, fp2, fp3, fp4, fp5, fp6, fp7;
fp1 = 1.5;
fp3 = ldexp( 1.0, -fp2.FractionBits );
fp2 = fp1 + fp3;
fp4 = fp3 / T( 2 );
fp5 = ldexp( 1.0, -fp2.FractionBits - 1 - bits );
fp6 = fp4 + fp5;
fp7 = fp1 + fp6;
cout << "Testing " << typeid( T::ValueType ).name() << " "
    << prec << "bit precision" << endl << endl;
cout    << "   " << fp1 << endl
    << " + " << fp3 << endl

    << " = " << fp2 << endl
    << endl;
cout    << "   " << fp3 << " / 2" << endl
    << " = " << fp4 << endl
    << endl;
cout    << "   " << fp4 << endl
    << " + " << fp5 << endl
    << " = " << fp6 << endl
    << endl;
cout    << "   " << fp1 << endl
    << " + " << fp6 << endl
    << " = " << fp7 << endl
    << ( fp7 == fp2 ? "OK." : "NG." ) << endl
    << endl;
cout << endl;
return fp7;
}
int main( void )
{
#ifdef __BORLANDC__
cout << "Compiled by Borland C++" << endl;
#endif
#ifdef _MSC_VER
cout << "Compiled by Visual C++" << endl;
#endif
SetNearRoundMode();
cout << "Testing nearest round mode" << endl << endl;
Test< Float >( Float::FractionBits, 24 );
Test< Float >( Float::FractionBits, 53 );
Test< Float >( Float::FractionBits, 64 );
Test< Double >( Double::FractionBits, 53 );
Test< Double >( Double::FractionBits, 64 );
Test< Double >( 63 - Double::FractionBits - 1, 64 );
Test< Double >( 63 - Double::FractionBits, 64 );
Test< LongDouble >( LongDouble::FractionBits, 64 );
return EXIT_SUCCESS;
}
///////////////////////////////////////////////////////////////////////////////

実行結果 7.1.4 のプログラムを PentiumIII 733MHz 上で実行した時の出力結果である。ただし、 “n bit precision” というのは内部の仮数部の計算精度である。また、出力形式は 2 進数であり、指数部の底も 2 である。

 Compiled by Borland C++
 Testing nearest round mode
 Testing float 24bit precision

  1.10000000000000000000000e+0
+ 1.00000000000000000000000e-23
= 1.10000000000000000000001e+0
  1.00000000000000000000000e-23 / 2
= 1.00000000000000000000000e-24
  1.00000000000000000000000e-24
+ 1.00000000000000000000000e-47
= 1.00000000000000000000001e-24
  1.10000000000000000000000e+0
+ 1.00000000000000000000001e-24
= 1.10000000000000000000001e+0
OK.

Testing float 53bit precision
  1.10000000000000000000000e+0
+ 1.00000000000000000000000e-23
= 1.10000000000000000000001e+0
  1.00000000000000000000000e-23 / 2
= 1.00000000000000000000000e-24
  1.00000000000000000000000e-24
+ 1.00000000000000000000000e-47
= 1.00000000000000000000001e-24
  1.10000000000000000000000e+0
+ 1.00000000000000000000001e-24
= 1.10000000000000000000001e+0
OK.

Testing float 64bit precision
  1.10000000000000000000000e+0
+ 1.00000000000000000000000e-23
= 1.10000000000000000000001e+0
  1.00000000000000000000000e-23 / 2
= 1.00000000000000000000000e-24
  1.00000000000000000000000e-24
+ 1.00000000000000000000000e-47
= 1.00000000000000000000001e-24
  1.10000000000000000000000e+0
+ 1.00000000000000000000001e-24
= 1.10000000000000000000001e+0
OK.

Testing double 53bit precision

  1.0000000000000000000000000000000000000000000000000000e-52 / 2
= 1.0000000000000000000000000000000000000000000000000000e-53
  1.0000000000000000000000000000000000000000000000000000e-53
+ 1.0000000000000000000000000000000000000000000000000000e-105
= 1.0000000000000000000000000000000000000000000000000001e-53
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000000000000000000000000000000000000000000000001e-53
= 1.1000000000000000000000000000000000000000000000000001e+0
OK.

Testing double 64bit precision
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000000000000000000000000000000000000000000000000e-52
= 1.1000000000000000000000000000000000000000000000000001e+0
  1.0000000000000000000000000000000000000000000000000000e-52 / 2
= 1.0000000000000000000000000000000000000000000000000000e-53
  1.0000000000000000000000000000000000000000000000000000e-53
+ 1.0000000000000000000000000000000000000000000000000000e-105
= 1.0000000000000000000000000000000000000000000000000001e-53
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000000000000000000000000000000000000000000000001e-53
= 1.1000000000000000000000000000000000000000000000000000e+0
NG.

Testing double 64bit precision
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000000000000000000000000000000000000000000000000e-52
= 1.1000000000000000000000000000000000000000000000000001e+0
  1.0000000000000000000000000000000000000000000000000000e-52 / 2
= 1.0000000000000000000000000000000000000000000000000000e-53
  1.0000000000000000000000000000000000000000000000000000e-53
+ 1.0000000000000000000000000000000000000000000000000000e-63
= 1.0000000001000000000000000000000000000000000000000000e-53
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000001000000000000000000000000000000000000000000e-53
= 1.1000000000000000000000000000000000000000000000000001e+0
OK.

Testing double 64bit precision
  1.1000000000000000000000000000000000000000000000000000e+0

+ 1.0000000000000000000000000000000000000000000000000000e-52
= 1.1000000000000000000000000000000000000000000000000001e+0
  1.0000000000000000000000000000000000000000000000000000e-52 / 2
= 1.0000000000000000000000000000000000000000000000000000e-53
  1.0000000000000000000000000000000000000000000000000000e-53
+ 1.0000000000000000000000000000000000000000000000000000e-64
= 1.0000000000100000000000000000000000000000000000000000e-53
  1.1000000000000000000000000000000000000000000000000000e+0
+ 1.0000000000100000000000000000000000000000000000000000e-53
= 1.1000000000000000000000000000000000000000000000000000e+0
NG.

Testing long double 64bit precision
  1.100000000000000000000000000000000000000000000000000000000000000e+0
+ 1.000000000000000000000000000000000000000000000000000000000000000e-63
= 1.100000000000000000000000000000000000000000000000000000000000001e+0
  1.000000000000000000000000000000000000000000000000000000000000000e-63 / 2
= 1.000000000000000000000000000000000000000000000000000000000000000e-64
  1.000000000000000000000000000000000000000000000000000000000000000e-64
+ 1.000000000000000000000000000000000000000000000000000000000000000e-127
= 1.000000000000000000000000000000000000000000000000000000000000001e-64
  1.100000000000000000000000000000000000000000000000000000000000000e+0
+ 1.000000000000000000000000000000000000000000000000000000000000001e-64
= 1.100000000000000000000000000000000000000000000000000000000000001e+0
OK.

7.2 斎藤卓也君の解答

私が課題のホームページを見たときには、 ’print the bits of double real number’ へのリンクが無く、自分で以下のプログラムを作り、 double の 64bit を書き出してみた。それぞれ、 near mode の浮動小数点演算で問題が出るという計算をさせ、その結果を調べてみた。

その後、そのホームページへの正しいリンクが貼られたので、実際にそこに載っていたプログラムでも実行してみたが、折角自分でプログラムを作ったので、ここではそれをメインに説明する。

7.2.1 作成したプログラム

/*

near mode の丸めが正しく行われるかどうかをチェックするプログラム

Programmed by Takuya Saito. 2001-4-17

*/

#include "stdio.h"
#include <float.h>
void Disp(volatile double x_n);
void FloatDisp(char* ptr);
void BitDisp(char c, int space_point);
int main()
{
unsigned int _controlftp(unsigned int new, unsigned int mask);
volatile double a,b,c,d;
/* near mode に設定 */
_controlfp(_RC_NEAR, _MCW_RC);
printf("a=");
a = 1.0/2/2/2/2/2/2/2/2/2/2
   /2/2/2/2/2/2/2/2/2/2
   /2/2/2/2/2/2/2/2/2/2
   /2/2/2/2/2/2/2/2/2/2
   /2/2/2/2/2/2/2/2/2/2
   /2/2/2;
Disp(a);
printf("b=");
b = a+1.0/2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2;
Disp(b);
printf("c=");

c = a+a/2/2/2/2/2/2/2/2/2/2
   /2;
Disp(c);
printf("d=");
d = c+1.0/2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2/2/2/2/2/2
     /2/2/2/2/2;
Disp(d);
printf("a+1=");
a=a+1.;
Disp(a);
printf("b+1=");
b=b+1.;
Disp(b);
printf("c+1=");
c=c+1.;
Disp(c);
printf("d+1=");
d=d+1.;
Disp(d);
printf("\n");
return 0;
}
/*

１０進数と IEEE754 の double を２進表示する関数

*/
void Disp(volatile double x_n)
{
    printf("%1.25g\n", x_n);
    FloatDisp((char*) &x_n);
    printf("\n");
}
/*

ptr で指定したポインタから、８バイト分を BitDisp 関数で２進数表示する関数

符号ビットと指数部、指数部と仮数部の間にはスペースを空けて表示する。

*/
void FloatDisp(char* ptr)
{
    int i;
    BitDisp(*(ptr+7),6);
    BitDisp(*(ptr+6),3);
    for (i=5; i>=0; i--){
        BitDisp(*(ptr+i),-1);
    }
    printf("\n");
}
/*

c で指定された１バイトの数字を２進表示する関数

space_point で、スペースを空けたい所を指定できる。

(-1 等を指定すると、スペースは空けない)

*/
void BitDisp(char c, int space_point)
{
    int i, wc;
    for(i=7;i>=0;i--) {
        wc = (c>>i) & 0x01;

if (i == space_point)
    printf(" ");
printf("%1d",wc);
}
}

7.2.2 実行結果

a=1.1102230246251565e-016
0 01111001010 0000000000000000000000000000000000000000000000000000
b=1.1102230246251568e-016
0 01111001010 0000000000000000000000000000000000000000000000000001
c=1.1107651257113993e-016
0 01111001010 0000000000100000000000000000000000000000000000000000
d=1.1107651257113995e-016
0 01111001010 0000000000100000000000000000000000000000000000000001
a+1=1
0 01111111111 0000000000000000000000000000000000000000000000000000
b+1=1.0000000000000002
0 01111111111 0000000000000000000000000000000000000000000000000001
c+1=1.0000000000000002
0 01111111111 0000000000000000000000000000000000000000000000000001
d+1=1.0000000000000002
0 01111111111 0000000000000000000000000000000000000000000000000001

Press any key to continue

２進数の表示は、符号ビット、指数部、仮数部がこの順でスペースで区切って表示するようにしてある。仮数部は正規化されているため、実際には先頭に 1 が一つ付いていると考える。

よって、例えば、 a に対して、 a + 1 の値は、

a + 1 = 1 + 1 × 2^{−53}

ということになり、これは、丁度 1 と 1 + 1 × 2^{−52} との中間の値であると考えることが出来る。そのため、 nearest mode で丸めた場合、 a + 1 の結果の様に、仮数部の最後のビット (LSB) が 0 であっても良い。

しかし、 b, c, d の値は、それぞれ a よりも僅かながら大きな値であるため、 b + 1, c + 1, d + 1 の値は、それぞれ 1 よりも 1 + 1 × 2^{−52} に近い値になる。

つまり、正常に nearest mode の丸めが機能するのであれば、 b+1, c+1, d+1 の値は、 1+1×2^{−52} とならなければならない。ということは、 b + 1, c + 1, d + 1 の結果の仮数部の最後のビット (LSB) が 1 でなければならない。

実際、上の結果を見ると、そのようになっているので、 nearest mode での浮動小数点演算は問題なく丸めが行われている。

課題のホームページにあった、 ’print the bits of double real number’ というリンク先の記事によると、 Pentium プロセッサには浮動小数点演算の nearest mode にバグがあり、正しい丸めが行われないとある。そして、そこにある実行結果を見ると、確かに b + 1 と c + 1 の計算結果の丸めが間違っている。

私は Intel Mobile celeron 600MHz 上で行ったので、この CPU では、浮動小数点演算の nearest mode の丸めには問題がないと言える。

そこで、さらに調べるために、現在私の手元にある他のコンピュータでも同様のプログラムを走らせ、試してみた。

実際に試した CPU は、 celeron 533MHz, Pentium 75MHz, Pentium 120MHz の３つである。だが、どの CPU でも、結果は Mobile celeron 600MHz 版と同じになり、 nearest mode の丸めは正常に行われた。試しに、私が作ったプログラムではなく、 ’print the bits of double real number’ のホームページに載っていた’rounding’ というプログラムでも試してみたが、残念ながら（？）どの CPU でもホームページに記載されているような丸めバグは発生しなかった。

以上の結果から、 ’print the bits of double real number’ というホームページで実際に計算させた Pentium がどのような物なのかが分からないが、初期の Pentium 75MHz や Pentium 120MHz、そして、最近の celeron 533MHz, Mobile celeron 600MHz には、 nearest mode にはバグが無いと言える。

バグがある CPU は、私が確認を行っていない MMX Pentium シリーズや、 Pentium II シリーズ等なのかも知れない。

7.3 匂坂岳志君の課題 3 の解答

例に挙げられたプログラムを実行して、 Intel Pentium III プロセッサと他のプロセッサの出力の違いをみた。

7.3.1 Intel Pentium III プロセッサについて

実験は第 5,6 端末室で行った。

                             Significand ([hidden bit]+52bit)
1+d= 1.00000000000000022e+00 [1.]0000 … 00000000 00000001

となっていることから、 1.10000000000000000e-16 を境界として、 IEEE754 の nearest mode （0 捨 1 入）による丸めを行えば、 1 + a, 1 + b, 1 + c, 1 + d それぞれの演算の結果は

                        Significand ([hidden bit]+52bit)
1.00000000000000000e+00 [1.]0000 … 00000000 00000000

または

                        Significand ([hidden bit]+52bit)
1.00000000000000022e+00 [1.]0000 … 00000000 00000001

となるはずである。

ここで、

          a= 1.11022302462515654e-16
          b= 1.11022302462515679e-16
          c= 1.11076512571139929e-16
          d= 1.11076512571139954e-16

となっていることから、結果は次のようになるはずである。

                             Significand ([hidden bit]+52bit)
1+a= 1.00000000000000022e+00 [1.]0000 … 00000000 0000000
1+b= 1.00000000000000022e+00 [1.]0000 … 00000000 00000001
1+c= 1.00000000000000022e+00 [1.]0000 … 00000000 00000001
1+d= 1.00000000000000022e+00 [1.]0000 … 00000000 00000001

しかし実際に実行してみると次のような結果が得られた。

                             Significand ([hidden bit]+52bit)
1+a= 1.00000000000000000e+00 [1.]0000 … 00000000 00000000
1+b= 1.00000000000000000e+00 [1.]0000 … 00000000 00000000
1+c= 1.00000000000000000e+00 [1.]0000 … 00000000 00000000
1+d= 1.00000000000000022e+00 [1.]0000 … 00000000 00000001

これは実験に使用した Pentium III プロセッサが、 IEEE754 の nearest mode （0 捨 1 入）の規定を厳密に満たしていないことを示している。

7.3.2 Intel 以外の他のプロセッサについて

Intel のプロセッサと同様に調べてみると、例に挙げられたプログラムでは 1+a, 1+b, 1+c, 1+d それぞれの演算のうち 1 + a のみ、規定を厳密に満たした演算結果でないようであった。

参考文献

[1] print the bits of double real number

   http://www.cmp.sanken.osaka-u.ac.jp/~nisimatu/computer
        /FloatingPointNumber/printbits.html

[2] 川口貴裕, IEEE Standard for Binary Floating-Point Arithmetic (和訳), 1998 http://www-amano.aa.cs.keio.ac.jp/members/kawaguti/memofpu/download/ieee754.ps

[3] GCC Home Page http://gcc.gnu.org/