浮動小数点数の指数部と仮数部への分解とその逆：logbとscalbn

浮動小数点数は指数部と仮数部の積で表される。

\[x=(-1)^s\cdot m\cdot 2^e\]

ここで、\(e\)は整数、\(m\)は\([1,b)\)または\([1/b,1)\)に正規化された数、あるいは\([b^{p-1},b^p)\)の範囲の整数である（流儀がいくつかあるということだ）。ただし\(b\)は基数で、\(p\)は仮数部の桁数とする。

プログラミング言語によっては、浮動小数点数を指数部と仮数部の組に変換したり、逆に指数部と仮数部の組から浮動小数点数を構築したりする関数を提供している。

ここではC言語を取り上げる。

C言語の規定 (C11)

C言語では ilogb/logb と scalbn という関数を提供している。

#include <math.h>
int ilogb(double x); // float版とlong double版もある
double logb(double x); // float版とlong double版もある
double scalbn(double x, int n); // float版とlong double版、それにnの型がlong intな版もある

ilogb と logb は指数部を返す。 scalbn は x を \(b^n\) 倍したものを返す。

このほか、frexp および ldexp という、「浮動小数点数を基数2と解釈したときの仮数部と指数部を返す」「\(2^n\) 倍する」関数もある。ほとんどの環境（FLT_RADIX == 2 の環境）では ldexp と scalbn は同一である。

#include <math.h>
double frexp(double value, int *exp); // float版とlong double版もある
double ldexp(double x, int exp); // float版とlong double版もある

ilogb / logb

#include <math.h>
int ilogb(double x); // float版とlong double版もある
double logb(double x); // float版とlong double版もある

ilogb と logb は指数部を返す。非正規化数は指数部の範囲を拡張して正規化したものとして扱う。

指数部\(e\)の流儀はいくつかあるが、ここでは\[\lvert x\rvert=m\cdot b^e\quad(1\leq m<b)\]とする。

ilogb の特殊な場合の対応は以下：

引数が0の場合、ilogb は FP_ILOGB0 という値を返す。これは INT_MIN または -INT_MAX のいずれかである。
引数が無限大の場合は INT_MAX を返す。
引数がNaNの場合は FP_ILOGBNAN という値を返す。これは INT_MIN または INT_MAX のいずれかである。
指数部を int で表現できない場合の返り値は未規定 (unspecified)。
Annex F準拠であれば、引数が0, 無限大、NaNの場合、あるいは正確な値を int で表現できない場合はinvalid例外が発生する。

logb の特殊な場合の対応は以下：

引数が0の場合、domain errorまたはpole errorが発生する。Annex F準拠であれば \(-\infty\) が返り、divide-by-zero例外が発生する。
引数が無限大の場合は \(+\infty\) が返る。
Annex F準拠であれば、（<math.h>の一般論として）引数がNaNの場合はNaNが返る。

scalbn

#include <math.h>
double scalbn(double x, int n); // float版とlong double版、それにnの型がlong intな版もある

scalbn は引数に与えられた浮動小数点数を \(b^n\) 倍する。

特殊な場合：

Annex F: \(x=\pm 0\)の場合は同じ符号の0が返ってくる。
Annex F: \(x=\pm\infty\)の場合は同じ符号の無限大が返ってくる。
Annex F:（<math.h>の一般論として）引数がNaNの場合はNaNが返る。

多くの場合、返ってくる値は正確である。しかし、返すべき値の指数部が小さすぎて正規化数で表現できない場合（アンダーフローする場合）は丸めが発生する可能性がある。

C11の段階ではアンダーフロー時の丸めについては規定がない。そのため、「数学的な値 \(x\cdot b^n\) に最も近い浮動小数点数」以外の値が返ってくる可能性がある。

実装上の注意点として、\(b^n\) を浮動小数点数として計算して掛け算する実装

double scalbn(double x, int n)
{
    double bn = pow(FLT_RADIX, n); // ここでは正確な値が得られると仮定する（本当は保証されていない）
    return x * bn;
}

はダメである。たとえば（doubleがbinary64な環境において） scalbn(0x1p-1074, 1074) は \(2^{-1074}\cdot 2^{1074}=1\) なので正確に1を返すべきであるが、\(2^{1074}\)は倍精度浮動小数点数で表現できないので上記の素朴な実装では無限大が返ってきてしまう。

少し改良して

double scalbn(double man, int exp)
{
    for (;;) {
        if (-1022 <= exp) {
            if (exp < 1024) {
                return man * pow(2.0, exp); // ここではpowで正確な値が得られると仮定する
            } else {
                man *= 0x1p1023;
                exp -= 1023;
            }
        } else {
            man *= 0x1p-1022;
            exp += 1022;
        }
    }
}

という実装がありうるが、これは結果が非正規化数の場合に「数学的な値 \(x\cdot b^n\) に最も近い浮動小数点数」以外の値が返ってくる可能性がある。まあC11の段階ではそれでもいいんだけど……。

具体例は x=0x1.555555555555bp-1, exp=-1024 である。数学的な値に最も近いのは 0x1.5555555555558p-1025 だが、上記の実装は 0x1.555555555556p-1025 を返す。

あなたが使っているlibc (libm)はどっちだろうか？筆者が使っているmacOSだと0x1.555555555556p-1025 が返ってきた。glibcやmuslだともうちょっと頑張っていると思う。

frexp

#include <math.h>
double frexp(double value, int *exp); // float版とlong double版もある

frexp は浮動小数点数を\([1/2,1)\)に正規化された仮数部と、基数2の指数部に分解する。

仮数部を正規化する範囲が違うので、（FLT_RADIX == 2 の場合）指数部の値は logb が返すものより1大きいことに注意しよう。

引数が特殊な値の場合は、

引数が0の場合は指数部も仮数部も0が返ってくる。
引数が非正規化数の場合も正規化が行われる。

さらにAnnex Fでは

引数が0の場合：同じ符号の0が返ってくる。
引数が無限大の場合：同じ符号の無限大が返ってくる。指数部の値は未規定 (unspecified) である。
引数がNaNの場合：何らかのNaNが返ってくる。指数部の値は未規定である。
C11の段階ではsignaling NaNに関する規定はない。

と規定している。

浮動小数点数の基数が2の場合、 frexp は丸めを起こさないことに注意しよう。返ってくる値は正確である。

ldexp

#include <math.h>
double ldexp(double x, int exp); // float版とlong double版もある

ldexp は浮動小数点数を \(2^{\mathtt{exp}}\) 倍したものを返す。注意点は scalbn と同様である。

IEEE 754の規定

IEEE 754-2019では浮動小数点数を指数部と仮数部に分解、あるいはその逆をするために scaleB および logB という演算が定義されている。

logB は

logBFormat logB(source)

という感じで、仮数部を \([1,b)\) に正規化した時の指数部を返す。logBFormat は浮動小数点数形式または整数形式である。

logBFormat が浮動小数点数の場合、

logB(NaN) はNaN
logB(∞) は \(+\infty\)
logB(0) は \(-\infty\) を返し、divideByZero例外を発生させる。

を返す。整数形式の場合は言語依存の値（絶対値が十分大きい）を返し、invalid operation例外を発生させる。

scaleB は

sourceFormat scaleB(source, logBFormat)

という型を持ち、 scaleB(x, N) は \(x\cdot b^N\) に最も近い浮動小数点数を返す。

IEEE 754-1985の時点ではscaleBに相当する演算は Scalb とされている（Annex A Recommended Functions and Predicates (Informative)）が、真の値を正確に表現できない場合の規定はない。

C言語とIEEE 754の関係

現行のC言語 (C11/C17) がAnnex Fで参照しているのはIEEE 754-1985相当なので、scalbn が真の値を正確に表現できない場合の規定がなくても仕方がないが、次期C標準（C23）のAnnex FではIEEE 754-2019を参照するようにアップデートされる。

そこではCの scalbn はIEEEのscaleB相当として言及されているので、真の値を正確に表現できない場合もscaleB相当の動作、つまり真の値に最も近い浮動小数点数を返すようにする必要があると思われる。

他の言語の場合

その辺のスクリプト言語だと指数部と仮数部への分解が利用できない場合がある。そういう場合は自分で実装する必要があるが、先に述べたように注意点がいくつかある。

指数部と仮数部への分解は、たとえばTypeScriptなら次のように書けるだろう：

function frexp(x: number): {man: number, exp: number}
{
    const a = Math.abs(x);
    if (a === 0.0 || a === Number.POSITIVE_INFINITY || Number.isNaN(a)) {
        return {man: x, exp: 0};
    } else {
        let e = Math.floor(Math.log2(a)) - 1; // 正確な値とは限らない
        for (;;) {
            const lower = 2**(e-1); // 正確な値が得られると仮定する（本当はimplementation-approximated valueでしかない）
            if (lower <= a) {
                if (a < lower * 2.0) { // lower * 2.0 は無限大かもしれない
                    return {man: x / lower * 0.5, exp: e};
                } else {
                    e += 1;
                }
            } else {
                e -= 1;
            }
        }
    }
}

逆は難しい。IEEEのscaleB準拠でなくて良いなら先に挙げたコードを移植すれば良いが、真の値に最も近い値を返したい場合はもうちょっと頑張る必要がある（どうやって？）。

Spread the love

雑記帳

人生やっていき

浮動小数点数の指数部と仮数部への分解とその逆：logbとscalbn

C言語の規定 (C11)

ilogb / logb

scalbn

frexp

ldexp

IEEE 754の規定

C言語とIEEE 754の関係

他の言語の場合

コメントを残すコメントをキャンセル

C言語の規定 (C11)

ilogb / logb

scalbn

frexp

ldexp

IEEE 754の規定

C言語とIEEE 754の関係

他の言語の場合

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル