nptclのブログ

Common Lisp処理系nptの開発メモです。https://github.com/nptcl/npt

非正規化浮動小数の有効桁数を求める

C言語

C言語にてCommon Lispのfloat-precision関数を作成したときのメモです。

非正規化の浮動小数とは、値が小さすぎるけど何とかして 浮動小数にて表現したいというギリギリの場合の方法です。

感覚としては、

1.23456 × 10の-103乗

みたいに表現したいけど指数-103は小さすぎて無理だから

0.00012 × 10の-99乗

のように表すのが非正規化です。

上記の例の場合、正規化の場合は1.23456という値から 6桁まで表現できているのがわかりますが、 非正規化の場合は0.000に押されて、12という2桁にまで減ってしまっています。 この「2桁」というのをC言語の世界で求める方法を説明します。

一応注意書きしておくと、上記の例は基数が10になっていますが、 実際はだいたい2なので注意してください。

まず普通の時、つまり正規化された浮動小数の桁数は、 C言語のfloat.hにdefineが用意されています。

FLT_MANT_DIG     24
DBL_MANT_DIG     53
LDBL_MANT_DIG    53, 64, 113とか？

もし浮動小数が正規化の場合は、何も考えずに上記の値を返却します。

では正規化か非正規化かを調べるにはどうしたらいいでしょうか。 思いつく方法は大小を比較することですが、 C99にはfpclassifyという関数が用意されています。 例えばdoubleを対象にした場合、

fpclassifyの戻り値が
  FP_NORMAL　なら　DBL_MANT_DIG

とします。 浮動小数点の値がゼロの場合はどうしましょうか？ ここでは有効桁数は0ということにします。

fpclassifyの戻り値が
  FP_ZEOR　なら　0

さて問題はfpclassifyの戻り値がFP_SUBNORMALの場合、つまり非正規化数の場合です。
返却値から先に示すと、

fpclassifyの戻り値が
  FP_SUBNORMAL　なら　DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(value);

となります。
DBL_MANT_DIGはすでに説明しました。
ilogb(value)は、指数部を整数で返却する関数です。

問題はDBL_MIN_EXPなのです。
DBL_MIN_EXPは何なのでしょうか。
検索すると2通りの意味が出てきます。

  • 基数のe乗で表現できる正規化数の最小値e
  • 基数のe-1乗で表現できる正規化数の最小値e

どっちなんだ。
つまり最小値そのものなのか、そうじゃなく最小値に1を足した数なのか？ 現物を調査すれば早いので結果を言いますと、e-1です。

IEEE754では、double型の正規化数の指数部の最小値は-1022になります。 つまり-1023に相当する値(binary値だと0になる)が非正規化数となります。

正規化数の最小値-1022より、

e-1 = -1022
e = -1021

となるので、

DBL_MIN_EXP = -1021

なのです。

では本当にC言語ではそうなのか？
コンパイラの問題だということもあり得ますので、 ちゃんと調べるためには規格書を見るしかありません。

ここであっているのかな。

Page 31
minimum negative integer such that FLT_RADIX raised to one less than that power is
 a normalized floating-point number, emin
FLT_MIN_EXP
DBL_MIN_EXP
LDBL_MIN_EXP

英語わからないです。 でもone less thanとは、「1を引いた数」という意味のようです。 なるほど、なぜ1を引いているのかは知りませんが、規格がそうなっているんですね。

ということなので、

FP_SUBNORMAL　なら　DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(value);

になります。

C言語にて、doubleのコードを示します。

#include <float.h>
#include <math.h>

int float_precision_double(double v)
{
    switch (fpclassify(v)) {
        case FP_NORMAL:
            return DBL_MANT_DIG;

        case FP_ZERO:
            return 0;

        case FP_SUBNORMAL:
            return DBL_MANT_DIG - (DBL_MIN_EXP - 1) + ilogb(v);

        default:
            fprintf(stderr, "error.\n");
            exit(1);
            return 0;
    }
}

実行の例を示します。

float_precision_double(123.456);
  →53
float_precision_double(0x01.0000p-1021);
  →53
float_precision_double(0x01.0000p-1022);
  →53
float_precision_double(0x01.0000p-1023);
  →52
float_precision_double(0x01.0000p-1024);
  →51
float_precision_double(0);
  →0

nptでは、ここで使っています。

最後に。
「非正規化」という名称はいろいろ種類があるように思えます。 デノーマル数という人もいれば、C言語では「SUBNORMAL」と表現されています。 Wikipediaの「非正規化数」が正しいのかもしれません。