💡
SIMD並列化ライブラリSmartVectorDotNet開発の知見まとめ(1) IEEE754浮動小数型の低レベル操作

2024/10/04に公開
現在、.Net用SIMD並列化ライブラリSmartVectorDotNetを開発しております。

本ライブラリの開発にあたっては数値型の低レベル操作、数学関数のアルゴリズムなど多くのノウハウを調査・実証し、パフォーマンスの追及を行ってきました。

ライブラリ自体はおおむね形になったので、後学のためその記録をここに残しておきます。

IEEE754浮動小数型の低レベル操作 (本稿)
SIMD演算の基礎
初等関数の実装
C#とdotnetの最適化

 binary32(float), binary64(double)の基礎知識現代的なコンピュータの多く、そして.Netの実行環境ではプリミティブな浮動小数型としてIEEE754のbinary32およびbinary64をサポートしています。
binary32の正規化数は次の数式で表すことができます。

\begin{array}{ll}
binary32 &= \mathrm{sign} \cdot (1+\mathrm{fraction}) \cdot 2^{\mathrm{exponent}-127}\\[1.6ex]
    &= (-1)^{b_{31}}
    \cdot \left( 1 + \sum_{i=0}^{22} b_i \cdot 2^{i-23} \right)
    \cdot 2^{-127 + \sum_{j=23}^{30}b_j 2^{j-23}}
\end{array}


(画像出典：IEEE754 - Wikipedia)
同様に、binary64の正規化数は次の数式で表すことができます。

\begin{array}{ll}
binary64 &= \mathrm{sign} \cdot (1+\mathrm{fraction}) \cdot 2^{\mathrm{exponent}-1023}\\[1.6ex]
    &= (-1)^{b_{63}}
    \cdot \left( 1 + \sum_{i=0}^{51} b_i \cdot 2^{i-52} \right)
    \cdot 2^{-1023 + \sum_{j=52}^{62}b_j 2^{j-52}}
\end{array}


(画像出典：IEEE754 - Wikipedia)
この変換を念頭に置くといくつかの数学関数を実装する際に効率化を行うことができます。例えばbinary32(b_{31}=0)の自然対数は

\begin{array}{ll}
\ln(binary32)
&= \ln \left\{ (1+\mathrm{fraction}) \cdot 2^{\mathrm{exponent}-127} \right\} \\[1.6ex]
&= \ln(1+\mathrm{fraction}) + \ln 2^{\mathrm{exponent}-127} \\[1.6ex]
&= \ln(1+\mathrm{fraction}) + (\mathrm{exponent}-127)\cdot\ln 2
\end{array}
\ln 2は定数になるので第2項は単純な四則演算のみで実装することができ、コア部分としては\ln(1+\mathrm{fraction})を実装すればよくなります。

後述しますが、0\le\mathrm{fraction}<1であることから\ln(1+\mathrm{fraction})はテイラー展開で計算できるため、正の正規化数全体に対して自然対数の実装を与えることができるようになります。

 剰余演算への応用binary32/binary64の性質をうまく利用すると剰余演算の効率化を図ることができます。
IEEE754において正の正規化数として表現可能な2つの実数

x = a \cdot 2^n (ただし1 \le a < 2, n\in \mathbb{Z})

y = b \cdot 2^m (ただし1 \le b < 2, m\in\mathbb{Z})
に対してx\ \mathrm{mod}\ yを考えます。
正の実数p, q, rに対して剰余演算\mathrm{mod}が

p\ \mathrm{mod}\ q = r \cdot \left(\frac{p}{r}\ \mathrm{mod}\ \frac{q}{r}\right)
の関係を満たすことを踏まえると、x\ \mathrm{mod}\ yは次のように変形できます。

n - m < 0ならば、x\ \mathrm{mod}\ y = x

n - m \ge 0ならば、x\ \mathrm{mod}\ y = 2^m\cdot(a\cdot 2^{n - m}\ \mathrm{mod}\ b)
ここで、IEEE754において

a\cdot 2^{n - m}はxの仮数部をマスクして整数型として再解釈し、さらにn - mだけ左シフトしたもの

（ただし、元のビット数に収まらない場合があるためオーバーフローを生じないよう十分にビット拡張されたものとする）

bはyの仮数部をマスクして整数型として再解釈したもの
つまり、オーバーフローが無視可能であるならば正規化浮動小数どうしの剰余演算も、整数の剰余演算とビットシフトに置換することができます。
また、仮数部をマスクしたということは、計算においては符号部・指数部のビットも利用することができます。このことを利用すればオーバーフロー対策として導入せざるを得なかったシフト→整数剰余計算の反復を最小限にとどめることができます。
以上により、IEEE754 binary32の剰余演算(正規化数限定)は以下のように実装できます。

このコードは入力に正規化数を要求する制限こそあるものの、C#標準の%演算子すら凌駕するほど高速に動作できます。
binary32正規化数剰余演算static float Modulo(float x, float y)
{
    const int signBitOffset = 31;
    const int expBitOffset = 23;
    const int signPartMask = 1 << signBitOffset;
    const int expPartMask = 0xFF << expBitOffset;
    const int fracPartMask = (1 << expBitOffset) - 1;
    const long economizedBit = 1L << expBitOffset;
    const int exponentBits = signBitOffset - expBitOffset;

    static void decompose(float x, out int sign, out int expo, out int frac)
    {
        var bin = Unsafe.As<float, int>(ref x);
        sign = (bin & signPartMask) >> signBitOffset;
        expo = (bin & expPartMask) >> expBitOffset;
        frac = bin & fracPartMask;
    }

    static float scale(int sign, int expo, int frac)
    {
        var bin = (sign << signBitOffset) | (expo << expBitOffset) | frac;
        return Unsafe.As<int, float>(ref bin);
    }

    decompose(x, out var s, out var n, out var a);
    decompose(y, out var _, out var m, out var b);
    var i = n - m;
    if (i < 0)
    {
        return x;
    }

    var c = (uint)(a | economizedBit);
    var d = (uint)(b | economizedBit);
    while (i > exponentBits)
    {
        c <<= exponentBits;
        c -= d * (c / d);
        i -= exponentBits;
    }
    c <<= (int)i;
    c -= d * (c / d);
    i = 0;

    if (c == 0)
    {
        return 0;
    }
    var eBitShift = BitOperations.LeadingZeroCount(c) - exponentBits;
    i -= eBitShift;
    c <<= eBitShift;
    return scale(s, m + i, (int)c & fracPartMask);
}
本番においては正規化数以外の処理もサポートする必要がありますが、このアルゴリズム自体は非常に有用であると考えられます。
binary32(float), binary64(double)の基礎知識

剰余演算への応用

Discussion
