備忘録: 浮動小数点データ

浮動小数点演算を実装するときに参照するパラメータや円周率など。

# 浮動小数点データ形式

数値 $N$ をいくつかの要素に分解して表現する方法

```math
\begin{eqnarray}
U & = & \left| N \right| \\
S & = & \frac{N}{U} \\
E & = & \lfloor log_B U \rfloor \\
T & = & \frac{U}{B^E} \\
N & = & S \times B^E \times T \\
\end{eqnarray} \\
```

で、$E$ と $T$ の範囲を定めることで、情報量を固定（32や64ビットなど）したまま幅広い数値を扱う。

| 項目 | 内容 | 備考 |
|:--|:--|:-:|
| 符号部 $S$ | 正なら $0$, 負なら $1$ | $sign\left(N\right)$ |
| 指数部 $E$ | 桁情報 | $floor \left( log \left( abs \left( N \right) \right) \div log \left( B \right) \right)$ |
| 仮数部 $T$ | 1.0 + 小数部 | $\left( abs\left(N\right) \div B^E \right)$ |

主に 2進形式 $\left(B=2\right)$ が使われているが、10進形式 $\left(B=10\right)$ も定義されている。

# ２進形式

各要素のビット並び

| | 符号部 $S$ | 　指数部 $E$　 | 　　　仮数部 $T$　　　 |
|:-:|:-:|:-:|:-:|
| ビット並び | $s$ | $ee \cdots ee$ | $\ \ mm \cdots mm$ |
| 正規数 ($E \neq 0$ の場合) | $\left(-1\right)^s$ | $2^{E-E_{bias}}$ | $1.mm \cdots mm$ |
| 非正規数 ($E = 0$ の場合) | $\left(-1\right)^s$ | $2^{-E_{bias}}$ | $0.mm \cdots mm$ |
| 無限大 | $\left(-1\right)^s$ | $11\cdots11$ | $0$ |
| 非数 $\left(NaN\right)$ | $\left(-1\right)^s$ | $11\cdots11$ | 非 $0$ |

データのビット幅と形式

| 形式 | $S$ | $E_{bits}$ | $T_{bits}$ | $E_{bias}$ |
|:-:|:-:|:-:|:-:|:-:|
| 単精度(32ビット) | 1 | 8 | 23 | 127 |
| 倍精度(64ビット) | 1 | 11 | 52 | 1023 |
| 拡張倍精度(80ビット)※ | 1 | 15 | 1+63 | 16383 |
| 四倍(?)精度(128ビット) | 1 | 15 | 112  | 16383 |

※拡張倍精度では、指数部 $E\neq0$ でも仮数部 $T$ に '1.0' 用の情報がある。

## 表現可能な精度

仮数部に入らない（切り捨てられる）値は、どんなものか。

| 精度 | 1.0 と次の数値の差分 |
|:-:|:--|
| $2^{-23}$ | 1.1920928955078125e-7 |
| $2^{-52}$ | 2.220446049250313080847263336181640625e-16 |
| $2^{-63}$ | 1.08420217248550443400745280086994171142578125e-19 |
| $2^{-112}$ | 1.925929944387235853055977942584927318538101648215388195239938795566558837890625e-34 |

<details><summary>1.0 と次の数値の実験プログラム </summary><div>

```c:sample.c
#include <stdio.h>

static float  fp32_1 = 1.0;
static float  fp32_2 = 1.0 + 1.1920928955078125e-7;
static double fp64_1 = 1.0;
static double fp64_2 = 1.0 + 2.220446049250313080847263336181640625e-16;

int main()
{
  printf("fp32_1: %#x\n", *(unsigned int *)(&fp32_1));
  printf("fp32_2: %#x\n", *(unsigned int *)(&fp32_2));
  printf("fp64_1: %#llx\n", *(unsigned long long *)(&fp64_1));
  printf("fp64_2: %#llx\n", *(unsigned long long *)(&fp64_2));
  return 0;
}
```

実行結果

```console:
$ ./sample
fp32_1: 0x3f800000
fp32_2: 0x3f800001
fp64_1: 0x3ff0000000000000
fp64_2: 0x3ff0000000000001
```

</div></details>

## 円周率π

円周率は、半分の $\pi \div 2$ が仮数部の形式に合っていて都合がよい値になります。

2倍や半分では、指数部が ±1 されます。

### π÷2 の2進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.10010010000111111011010 <br/> 1.10010010000111111011011 |
| 倍精度 | 1.1001001000011111101101010100010001000010110100011000 <br/> 1.1001001000011111101101010100010001000010110100011001 |
| 拡張倍精度 | 1.100100100001111110110101010001000100001011010001100001000110100 <br/> 1.100100100001111110110101010001000100001011010001100001000110101 |
| 四倍精度 | 1.1001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111000 <br/> 1.1001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111001 |

### π÷2 の10進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.5707962512969970703125 <br/> 1.57079637050628662109375 |
| 倍精度 | 1.5707963267948965579989817342720925807952880859375 <br/> 1.5707963267948967800435866593034006655216217041015625 |
| 拡張倍精度 | 1.5707963267948966191479842624545426588156260550022125244140625 <br/> 1.570796326794896619256404479703093102216371335089206695556640625 |
| 四倍精度 | 1.5707963267948966192313216916397513987395340490686477865022521659371483593314877680313657037913799285888671875 <br/> 1.5707963267948966192313216916397515913325284877922330921000464244298802131416525895701852277852594852447509765625 |

### π÷2 の16進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.921fb4 <br/> 1.921fb6 |
| 倍精度 | 1.921fb54442d18 <br/> 1.921fb54442d19 |
| 拡張倍精度 | 1.921fb54442d18468 <br/> 1.921fb54442d18469 |
| 四倍精度 | 1.921fb54442d18469898cc51701b8 <br/> 1.921fb54442d18469898cc51701b9 |

## ネイピア数e

ネイピア数も、半分の $e \div 2$ が仮数部の形式に合っていて都合がよい値になります。

2倍や半分では、指数部が ±1 されます。

### e÷2 の2進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.01011011111100001010100 <br/> 1.01011011111100001010101 |
| 倍精度 | 1.0101101111110000101010001011000101000101011101101001 <br> 1.0101101111110000101010001011000101000101011101101010 |
| 拡張倍精度 | 1.010110111111000010101000101100010100010101110110100101010011010 <br/> 1.010110111111000010101000101100010100010101110110100101010011011 |
| 四倍精度 | 1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111010 <br/> 1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111011 |

### e÷2 の10進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.359140872955322265625 <br/> 1.35914099216461181640625 |
| 倍精度 | 1.3591409142295225453977991492138244211673736572265625 <br> 1.359140914229522767442404074245132505893707275390625 |
| 拡張倍精度 | 1.35914091422952261760566383674841972606373019516468048095703125 <br/> 1.359140914229522617714084053996970169464475475251674652099609375 |
| 四倍精度 | 1.359140914229522617680143735676331157179210933596774431334615430163833584009669408487752662040293216705322265625 <br/> 1.3591409142295226176801437356763313497722053723203597369324096886565654378198342300265721860341727733612060546875 |

### e÷2 の16進数

| 精度 | 切り下げ値 <br/> 切り上げ値 |
|:-:|:--|
| 単精度 | 1.5bf0a8 <br/> 1.5bf0aa |
| 倍精度 | 1.5bf0a8b145769 <br> 1.5bf0a8b14576a |
| 拡張倍精度 | 1.5bf0a8b145769534 <br/> 1.5bf0a8b145769536 |
| 四倍精度 | 1.5bf0a8b1457695355fb8ac404e7a <br/> 1.5bf0a8b1457695355fb8ac404e7b |

# 10進形式

10進形式では、指数部に仮数部の一部を混ぜて符号化した複合部がある。さらに、小数部は10進数3桁を10ビットで符号化される。

### 各要素のビット並び

| | 符号部 $S$ | 　　　　　　複合部 $G$　　　　　　 <br/> 指数部 $E$ と先頭桁 | 仮数部 $T$ |
|:-:|:-:|:-:|:-:|
| 　　ビット並び　　 | $s$ | $ww \cdots ww$ | $tt \cdots tt$ |
| ※ 先頭桁 <br/> $0\cdots7=ddd$ | $\left(-1\right)^s$ | $e_0e_1ddd\ ee \cdots ee$ <br/> $e_0e_1 \neq 11$ | 小数部 |
| 先頭桁 <br/> $8,9=8+d$ | $\left(-1\right)^s$ | $11e_0e_1d\ ee \cdots ee$ <br/> $e_0e_1 \neq 11$ | 小数部 |
| 無限大 | $\left(-1\right)^s$ | $11110\ xxx \cdots xx$ | 無効 |
| 非数 $\left(qNaN\right)$ | $\left(-1\right)^s$ | $11111\ 0xx \cdots xx$ | 識別用 |
| 非数 $\left(sNaN\right)$ | $\left(-1\right)^s$ | $11111\ 1xx \cdots xx$ | 識別用 |

※仮数部が "$0.00\cdots$" ならば指数部に関係なくゼロになる。

### データのビット幅と形式

| decimal <br/> $k$ | $S$ <br/> $1$ |　　 $G_{bits}$　　 <br/> $5+w$ | $T_{bits}$ <br/> $t$ | 　桁数 $p$　 <br/> &nbsp; | $E_{bias}$ <br/> &nbsp; |
|:-:|:-:|:-:|:-:|:-:|:-:|
| $32$ | $1$ | $5+6$ | $20$ | $1+6$ | $101$ |
| $64$ | $1$ | $5+8$ | $50$ | $1+15$ | $398$ |
| $96$ | $1$ | $5+10$ | $80$ | $1+24$ | $1559$ |
| $128$ | $1$ | $5+12$ | $110$ | $1+33$ | $6176$ |

指数部 $E = e_0e_1\ ee \cdots ee $ の値は $e_0e_1 \neq 11$ なので範囲は
　　$0$ 〜 $\left( 3 \times 2^w - 1 \right)$
になる。

### 小数部の符号化

10進数3桁の10ビット符号化形式

| 　　　　　　形式　　　　　　 | 　$d_1\times100$　 | 　$d_2\times10$　 | 　$d_3\times1$　 |
|:-:|:-:|:-:|:-:|
| $[b_0 b_1 b_2][b_3 b_4 b_5]\ 0\ [b_7 b_8 b_9]$ <br/> $[d_1 d_1 d_1][d_2 d_2 d_2]\ 0\ [d_3 d_3 d_3]$ | $b_0 b_1 b_2$ <br/> $0\cdots7$ | $b_3 b_4 b_5$ <br/> $0\cdots7$ | $b_7 b_8 b_9$ <br/> $0\cdots7$ |
| $[b_0 b_1 b_2][b_3 b_4 b_5]\ 100\ [b_9]$ <br/> $[d_1 d_1 d_1][d_2 d_2 d_2]\ 100\ [d_3]$ | $b_0 b_1 b_2$ <br/> $0\cdots7$ | $b_3 b_4 b_5$ <br/> $0\cdots7$ | $100b_9$ <br/> $8+b_9$ |
| $[b_0 b_1 b_2][b_3 b_4][b_5]\ 101\ [b_9]$ <br/> $[d_1 d_1 d_1][d_3 d_3][d_2]\ 101\ [d_3]$ | $b_0 b_1 b_2$ <br/> $0\cdots7$ | $100b_5$ <br/> $8+b_5$ | $b_3 b_4 b_9$ <br/> $0\cdots7$  |
| $[b_0 b_1][b_2][b_3 b_4 b_5]\ 110\ [b_9]$ <br/> $[d_3 d_3][d_1][d_2 d_2 d_2]\ 110\ [d_3]$ | $100b_2$ <br/> $8+b_2$ | $b_3 b_4 b_5$ <br/> $0\cdots7$ | $b_0 b_1 b_9$ <br/> $0\cdots7$ |
| $[b_0 b_1][b_2]\ 00\ [b_5]\ 111\ [b_9]$ <br/> $[d_3 d_3][d_1]\ 00\ [d_2]\ 111\ [d_3]$ | $100b_2$ <br/> $8+b_2$ | $100b_5$ <br/> $8+b_5$ | $b_0 b_1 b_9$ <br/> $0\cdots7$ |
| $[b_0 b_1][b_2]\ 01\ [b_5]\ 111\ [b_9]$ <br/> $[d_2 d_2][d_1]\ 01\ [d_2]\ 111\ [d_3]$ | $100b_2$ <br/> $8+b_2$ | $b_0 b_1 b_5$ <br/> $0\cdots7$ | $100b_9$ <br/> $8+b_9$ |
| $[b_0 b_1 b_2]\ 10\ [b_5]\ 111\ [b_9]$ <br/> $[d_1 d_1 d_1]\ 10\ [d_2]\ 111\ [d_3]$ | $b_0 b_1 b_2$ <br/> $0\cdots7$ | $100b_5$ <br/> $8+b_5$ | $100b_9$ <br/> $8+b_9$ |
| $xx\ [b_2]\ 11\ [b_5]\ 111\ [b_9]$ <br/> $xx\ [d_1]\ 11\ [d_2]\ 111\ [d_3]$ | $100b_2$ <br/> $8+b_2$ | $100b_5$ <br/> $8+b_5$ | $100b_9$ <br/> $8+b_9$ |

<details><summary>符号化／復号化 プログラム</summary><div>

```c:packed_decimal.c
#include <stdio.h>

unsigned int decode_packed_decimal(unsigned int e)
{
  unsigned int d1 = (e >> 7) & 7;
  unsigned int d2 = (e >> 4) & 7;
  unsigned int d3 = (e >> 0) & 7;
  unsigned int m1 = d1 & 6;
  unsigned int m2 = d2 & 6;
  unsigned int m3 = d3 & 6;
  unsigned int l1 = d1 & 1;
  unsigned int l2 = d2 & 1;
  unsigned int l3 = d3 & 1;

  switch ((e >> 1) & 7)
    {
    default: break;
    case 0x04:
      d3 = l3 | 8;
      break;

    case 0x05:
      d2 = l2 | 8;
      d3 = l3 | m2;
      break;

    case 0x06:
      d1 = l1 | 8;
      d3 = l3 | m1;
      break;

    case 0x07:
      switch ((e >> 5) & 3)
	{
	case 0:
	  d1 = l1 | 8;
	  d2 = l2 | 8;
	  d3 = l3 | m1;
	  break;

	case 1:
	  d1 = l1 | 8;
	  d2 = l2 | m1;
	  d3 = l3 | 8;
	  break;

	case 2:
	  d2 = l2 | 8;
	  d3 = l3 | 8;
	  break;

	case 3:
	  d1 = l1 | 8;
	  d2 = l2 | 8;
	  d3 = l3 | 8;
	  break;
	}
    }
  return (d1 * 100 + d2 * 10 + d3);
}

unsigned int encode_packed_decimal(unsigned int n)
{
  unsigned int d1 = (n / 100) % 10;
  unsigned int d2 = (n /  10) % 10;
  unsigned int d3 = (n /   1) % 10;

  unsigned int h1 = (d1 & 8);
  unsigned int h2 = (d2 & 8);
  unsigned int h3 = (d3 & 8);
  /* unsigned int m1 = d1 & 6; */
  unsigned int m2 = d2 & 6;
  unsigned int m3 = d3 & 6;
  unsigned int l1 = d1 & 1;
  unsigned int l2 = d2 & 1;
  unsigned int l3 = d3 & 1;

  d1 <<= 7; l1 <<= 7;
  d2 <<= 4; l2 <<= 4;
  /* d3 <<= 0; l3 <<= 0; */

  switch ((h1 >> 1) | (h2 >> 2) | (h3 >> 3))
    {
    case 0: return 0x00 | d1 | d2 | d3;
    case 1: return 0x08 | d1 | d2 | l3;
    case 2: return 0x0a | d1 | l2 | l3 | (m3 << 4);
    case 4: return 0x0c | l1 | d2 | l3 | (m3 << 7);
    case 6: return 0x0e | l1 | l2 | l3 | (m3 << 7);
    case 5: return 0x2e | l1 | l2 | l3 | (m2 << 7);
    case 3: return 0x4e | d1 | l2 | l3;
    case 7: return 0x6e | l1 | l2 | l3;
    }
  return -1;
}

int main(int argc, char **argv)
{
  unsigned int n;

  for (n = 0; n < 1000; n++)
    {
      unsigned int e = encode_packed_decimal(n);
      unsigned int d = decode_packed_decimal(e);
      printf("%3d > %03x > %3d - %s\n", n, e, d, (n == d) ? "ok" : "ng");
    }
  return 0;
}
```

</div></details>

# 2進数と10進

有限の小数点以下 $b$ 桁の 2進数では、10進数の "0.1" は表現できない

```math
\frac{A}{2^b} \neq 0.1
```

ので、小数点以下 $d$ 桁の正確な 10進数演算が欲しいときがあります。

有限の小数点以下 $b$ 桁の 2進数は、有限の小数点以下 $b$ 桁の 10進数で

```math
\begin{eqnarray} \\
\frac{A}{\left( 2 \times 5 \right)^b} & = & \frac{C}{2^b} \\
A & = & 5^b C \\
\end{eqnarray}
```

誤差のない表現にできますが、2進数と同じ桁数表示になります。

| 基数 | 値 |
|:-:|:-:|
| 2進数 | $1.10010010000111111011011$ |
| 10進数 | $1.57079637050628662109375$  |

プログラムを組んでいると、10進数の面倒臭さを感じることが多い…

---