浮動小数点数について復習してみた

最近、当たり前のことを当たり前のことのように理解していることが、とても大切なことのように感じています。その中でも浮動小数点数は日ごろ当たり前のように触っているのに、実際どういう風に動いているのか全く理解できていないことに気づきびっくりしました。そこで今回は浮動小数点数について復習してみました。

# 浮動小数点数とは

浮動小数点数はコンピュータで `3.14` のような実数を近似して表現する方法です。かつてはいろいろな浮動小数点数の表現方法があったようなのですが、我々が日ごろ扱う環境ではたいていの場合 IEEE 754 とよばれる国際規格に基づいた表現方法が採用されています。

私は最近は JavaScript を扱うことが多いのですが、JavaScript では数値 (`Number`) は IEEE 754 の倍精度浮動小数点数というフォーマットで表現されています。倍精度浮動小数点数では値を符号部1ビット, 指数部11ビット、仮数部52ビットの合計64ビットで表現します。

JavaScript では [`ArrayBuffer`](https://developer.mozilla.org/ja/docs/Web/JavaScript/Reference/Global_Objects/ArrayBuffer) と [`DataView`](https://developer.mozilla.org/ja/docs/Web/JavaScript/Reference/Global_Objects/DataView) を利用すると値のバイト列を取り出すことができます。以下は数値 `3.0` を16進表記で表示した例です。

```js
function numberToHexString(value) {
    // 8バイトの領域を作る
    const buffer = new ArrayBuffer(8);

    // 0バイト目に引数の値を倍精度浮動小数点数としてリトルエンディアンで書き込む。
    const dataView = new DataView(buffer);
    dataView.setFloat64(0, value, true);

    // 16進数文字列に変換する。
    let result = "";
    for (let i = 0; i < 8; i++) {
        const uint8 = dataView.getUint8(i);
        if (uint8 <= 0x0F) result += "0"
        result += uint8.toString(16) + " ";
    }
    return result;
}

console.log("+0 +1 +2 +3 +4 +5 +6 +7")
console.log(numberToHexString(3.0));
```

実行するとコンソールに次のように表示されます。

```shell-session
+0 +1 +2 +3 +4 +5 +6 +7
00 00 00 00 00 00 08 40
```

リトルエンディアンで格納したのでアドレスが大きいほうが上位バイトです。このままではよくわからないので、符号部、指数部、仮数部に分解すると次のようになります。

- 符号部 (1ビット): `0x00 (0b0)`
- 指数部 (11ビット): `0x400 (0b10000000000)`
- 仮数部 (52ビット): `0x8000000000000 (0b1000..(0が48個))`

符号部の 0 は値が正の数であることを表しています。指数部の `0x400` は2を基数としたときの指数の値に1023を足した値が入っています。`0x400 - 1023` より指数の値は `1` です。仮数部は整数部が1ビットの固定小数点数 (Q1) の小数部が入っています。つまり、2進数の小数 `0b1.1` 10進数だと `1.5` が入っています。 仮数の値を2の指数乗すると `3.0` になり予想通りの値になっていることがわかります。

# 浮動小数点数の乗算

これでビット列で数値を表現できるようになりました。次は浮動小数点数同士の乗算について考えます。例えば `3.0` と `4.0` を掛けて `12.0` になることをビット列で考えます。

`3.0` と `4.0` のビット列はそれぞれ次のようになります。

|      |16進表現             |符号部|指数部  |仮数部           |
|------|--------------------|------|-------|-----------------|
|`3.0` |`0x4008000000000000`|`0`   |`0x400`|`0x8000000000000`|
|`4.0` |`0x4010000000000000`|`0`   |`0x401`|`0x0000000000000`|
|`12.0`|`0x4028000000000000`|`0`   |`0x402`|`0x8000000000000`|

ここでは `x = 3.0, y = 4.0` として、`x * y` を計算します。まずは結果の符号を考えます。これは右辺の符号と左辺の符号が等しければ結果は正の数、異なっていれば負の数です。`!=` 演算子か `^` (XOR) 演算子で計算できます。

今回は両辺どちらも `0` と等しいので結果の符号部は `0` になります。

次に指数部を考えます。学校の数学の時間で教わったように、2つの値を掛け算するとき、指数の値は足し合わせればよいので足し合わせます。倍精度浮動小数点数の指数部は 1023 加えられているので、

```
0x400 - 1023 + 0x401 - 1023 + 1023 = 0x402
```

より結果の指数部は `0x402` になります。

最後に仮数部です。整数部が1ビット、小数部が52ビットの値同士を掛け算すると結果は、整数部2ビット小数部104ビットの値になります。

```
0x18000000000000 * 0x10000000000000 = 0x180000000000000000000000000
```

結果の整数部は `0x1` です。ここで2ビット目が1の場合は指数部と仮数部を調整する必要がありますが、今回は不要です。小数部の上から52ビットを取り出して、結果は `0x8000000000000` になります。

これまで計算してきた符号部と指数部と仮数部を連結すると、結果の倍精度浮動小数点数表現が得られます。

# 思い付き

実際に処理を追っていくと、浮動小数点数の乗算は整数の加算(指数部)と整数の乗算(仮数部)の組み合わせであることがわかります。ということは、適切に値の格納と抽出ができれば、浮動小数点数の乗算で、加算と乗算を一度に行えることになります。

ということで、C言語で浮動小数点数演算を使って整数の加算と乗算を同時に処理するプログラムを書いてみました。

まずは以下の構造体を作ります。

```c
typedef union {
    uint64_t v;
    double f;
} f64;
```

2つの整数 `v1` と `v2` を一つの倍精度浮動小数点数に埋め込みます。`v1` を加算 `v2` を乗算に使います。構造上の制限から、`v1` は9ビット、`v2` は13ビットまでの符号なし整数値を埋め込むことができます。2つの値を倍精度浮動小数点数に埋め込むには次のようにします。

```c
// v1: 9ビット, v2: 13ビット
double pack(uint16_t v1, uint16_t v2)
{
    return (f64){ .v = ((v1 + 1023ul) << 52ul) | (v2 << 26ul) }.f; 
}
```

それぞれ値を取り出す操作は次のようになります。

```c
uint16_t unpack1(double f)
{
    return (((f64){ .f = f }.v >> 52) - 1023) & 0x1FF;
}
```

```c
uint16_t unpack2(double f)
{
    return (f64){ .f = f }.v & 0x1FFF;
}
```

以下のプログラムを実行すると、画面に `5 6` と表示されます。つまり `f1 * f2` で `2 + 3` と `2 * 3` が同時に計算できたことがわかります。

```c
int main(void)
{    
    double f1 = pack(2, 2);
    double f2 = pack(3, 3);

    double f3 = f1 * f2;
    
    printf("%d %d\n", unpack1(f3), unpack2(f3)); // => 5 6
    
    return 0;
}
```

ちなみに、埋め込み方法を変えると、`v2` で行われる処理を足し算に変更できます。足し算の場合は最大26ビットまで扱えます。

```
double pack2(uint16_t v1, uint16_t v2)
{
    return (f64){ .v = ((v1 + 1023ul) << 52ul) | v2 }.f; 
}
```

# プログラム例

以下は `n` の階乗を計算するプログラム例です。階乗を計算するための乗算の処理と、ループカウンタのインクリメントをひとつの乗算演算子で処理している点がポイントです。

```c
int fact(int n)
{
    const double c = 1 + DBL_EPSILON;
    double p1 = c, p2 = c * 2;

    while (p1 < pack(n - 1, 0)) {
        p2 *= c;
        p1 = shift(p1) * shift(p2);
    }

    return unpack2(p1);
}
```

ここで注意しなければならないのは、`pack()` で作成した値同士を2つかけると `pack2()` 相当の値になる点である。そこで `pack2` 値を `pack1` 値に変換する関数 `shift` を作成しました。

```c
double shift(double f)
{
    uint64_t v = (f64){ .f = f }.v;
    return (f64){ .v = (v & (0xFFF0ul << 48)) | ((v & 0xFFFF) << 26) }.f;
}
```

処理の内容は、単に `v2` 側の値を左に26ビットシフトしているだけです。

main 関数を作って実行したら、以下のように表示されました。ちゃんと計算できているようです。

```
fact(0) = 1
fact(1) = 1
fact(2) = 2
fact(3) = 6
fact(4) = 24
fact(5) = 120
fact(6) = 720
fact(7) = 5040
```

# まとめ

なんか話がどんどん脱線していきましたが、浮動小数点数について復習していたら、浮動小数点数の乗算で整数の加算と乗算を一気に処理できることに気づきました。特にこの方法自体は何の役にも立ちそうにありませんが、浮動小数点数の演算がなんとなく身近になった気がします。これからは自信を持って JavaScript で数値を扱えそうです。

# 参考文献
- [浮動小数点数 (Wikipedia)](https://ja.wikipedia.org/wiki/%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%95%B0)
- [IEEE 754 (Wikipedia)](https://ja.wikipedia.org/wiki/IEEE_754)