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

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

浮動小数点データ形式

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

\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

1.0 と次の数値の実験プログラム

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;
}

実行結果

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

円周率π

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

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

π÷2 の2進数

精度	切り下げ値 切り上げ値
単精度	1.10010010000111111011010 1.10010010000111111011011
倍精度	1.1001001000011111101101010100010001000010110100011000 1.1001001000011111101101010100010001000010110100011001
拡張倍精度	1.100100100001111110110101010001000100001011010001100001000110100 1.100100100001111110110101010001000100001011010001100001000110101
四倍精度	1.1001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111000 1.1001001000011111101101010100010001000010110100011000010001101001100010011000110011000101000101110000000110111001

π÷2 の10進数

精度	切り下げ値 切り上げ値
単精度	1.5707962512969970703125 1.57079637050628662109375
倍精度	1.5707963267948965579989817342720925807952880859375 1.5707963267948967800435866593034006655216217041015625
拡張倍精度	1.5707963267948966191479842624545426588156260550022125244140625 1.570796326794896619256404479703093102216371335089206695556640625
四倍精度	1.5707963267948966192313216916397513987395340490686477865022521659371483593314877680313657037913799285888671875 1.5707963267948966192313216916397515913325284877922330921000464244298802131416525895701852277852594852447509765625

π÷2 の16進数

精度	切り下げ値 切り上げ値
単精度	1.921fb4 1.921fb6
倍精度	1.921fb54442d18 1.921fb54442d19
拡張倍精度	1.921fb54442d18468 1.921fb54442d18469
四倍精度	1.921fb54442d18469898cc51701b8 1.921fb54442d18469898cc51701b9

ネイピア数e

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

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

e÷2 の2進数

精度	切り下げ値 切り上げ値
単精度	1.01011011111100001010100 1.01011011111100001010101
倍精度	1.0101101111110000101010001011000101000101011101101001 1.0101101111110000101010001011000101000101011101101010
拡張倍精度	1.010110111111000010101000101100010100010101110110100101010011010 1.010110111111000010101000101100010100010101110110100101010011011
四倍精度	1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111010 1.0101101111110000101010001011000101000101011101101001010100110101010111111011100010101100010000000100111001111011

e÷2 の10進数

精度	切り下げ値 切り上げ値
単精度	1.359140872955322265625 1.35914099216461181640625
倍精度	1.3591409142295225453977991492138244211673736572265625 1.359140914229522767442404074245132505893707275390625
拡張倍精度	1.35914091422952261760566383674841972606373019516468048095703125 1.359140914229522617714084053996970169464475475251674652099609375
四倍精度	1.359140914229522617680143735676331157179210933596774431334615430163833584009669408487752662040293216705322265625 1.3591409142295226176801437356763313497722053723203597369324096886565654378198342300265721860341727733612060546875

e÷2 の16進数

精度	切り下げ値 切り上げ値
単精度	1.5bf0a8 1.5bf0aa
倍精度	1.5bf0a8b145769 1.5bf0a8b14576a
拡張倍精度	1.5bf0a8b145769534 1.5bf0a8b145769536
四倍精度	1.5bf0a8b1457695355fb8ac404e7a 1.5bf0a8b1457695355fb8ac404e7b

10進形式

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

各要素のビット並び

	符号部 $S$	複合部 $G$　　　　　　 指数部 $E$ と先頭桁	仮数部 $T$
ビット並び	$s$	$ww \cdots ww$	$tt \cdots tt$
※ 先頭桁 $0\cdots7=ddd$	$\left(-1\right)^s$	$e_0e_1ddd\ ee \cdots ee$ $e_0e_1 \neq 11$	小数部
先頭桁 $8,9=8+d$	$\left(-1\right)^s$	$11e_0e_1d\ ee \cdots ee$ $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 $k$	$S$ $1$	$G_{bits}$　　 $5+w$	$T_{bits}$ $t$	桁数 $p$	$E_{bias}$
$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]$ $[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$ $0\cdots7$	$b_3 b_4 b_5$ $0\cdots7$	$b_7 b_8 b_9$ $0\cdots7$
$[b_0 b_1 b_2][b_3 b_4 b_5]\ 100\ [b_9]$ $[d_1 d_1 d_1][d_2 d_2 d_2]\ 100\ [d_3]$	$b_0 b_1 b_2$ $0\cdots7$	$b_3 b_4 b_5$ $0\cdots7$	$100b_9$ $8+b_9$
$[b_0 b_1 b_2][b_3 b_4][b_5]\ 101\ [b_9]$ $[d_1 d_1 d_1][d_3 d_3][d_2]\ 101\ [d_3]$	$b_0 b_1 b_2$ $0\cdots7$	$100b_5$ $8+b_5$	$b_3 b_4 b_9$ $0\cdots7$
$[b_0 b_1][b_2][b_3 b_4 b_5]\ 110\ [b_9]$ $[d_3 d_3][d_1][d_2 d_2 d_2]\ 110\ [d_3]$	$100b_2$ $8+b_2$	$b_3 b_4 b_5$ $0\cdots7$	$b_0 b_1 b_9$ $0\cdots7$
$[b_0 b_1][b_2]\ 00\ [b_5]\ 111\ [b_9]$ $[d_3 d_3][d_1]\ 00\ [d_2]\ 111\ [d_3]$	$100b_2$ $8+b_2$	$100b_5$ $8+b_5$	$b_0 b_1 b_9$ $0\cdots7$
$[b_0 b_1][b_2]\ 01\ [b_5]\ 111\ [b_9]$ $[d_2 d_2][d_1]\ 01\ [d_2]\ 111\ [d_3]$	$100b_2$ $8+b_2$	$b_0 b_1 b_5$ $0\cdots7$	$100b_9$ $8+b_9$
$[b_0 b_1 b_2]\ 10\ [b_5]\ 111\ [b_9]$ $[d_1 d_1 d_1]\ 10\ [d_2]\ 111\ [d_3]$	$b_0 b_1 b_2$ $0\cdots7$	$100b_5$ $8+b_5$	$100b_9$ $8+b_9$
$xx\ [b_2]\ 11\ [b_5]\ 111\ [b_9]$ $xx\ [d_1]\ 11\ [d_2]\ 111\ [d_3]$	$100b_2$ $8+b_2$	$100b_5$ $8+b_5$	$100b_9$ $8+b_9$

符号化／復号化 プログラム

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;
}

2進数と10進

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

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

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

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

\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進数の面倒臭さを感じることが多い…

---