浮動小数点データから短い10進数表記を求める

コンピューターで小数点以下を扱う数値計算は、CPU に演算機能が搭載されている [IEEE 754 浮動小数点](https://qiita.com/ikiuo/items/6c02524a35f676ec753b)がよく使われています。ここでは様々な開発言語で使用されている IEEE 754 倍精度浮動小数点型(F64と略記)のうち、正の数を前提とします。

算出した F64 データの全精度で表記するとき、16進浮動小数点表記ができる言語は長くても20数文字です。16進表記できない言語では10進数に頼ることになります。すると、2進数での小数点以下の桁数が10進数の桁数になるので小さい値ほど長くなります。非正規化数の正の最小値 $2^{-1074}$ を厳密値にすると下1074桁になります。

```text:全桁(下1074桁)
0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004940656458412465441765687928682213723650598026143247644255856825006755072702087518652998363616359923797965646954457177309266567103559397963987747960107818781263007131903114045278458171678489821036887186360569987307230500063874091535649843873124733972731696151400317153853980741262385655911710266585566867681870395603106249319452715914924553293054565444011274801297099995419319894090804165633245247571478690147267801593552386115501348035264934720193790268107107491703332226844753335720832431936092382893458368060106011506169809753078342277318329247904982524730776375927247874656084778203734469699533647017972677717585125660551199131504891101451037862738167250955837389733598993664809941164205702637090279242767544565229087538682506419718265533447265625
```

```text:指数表現形式(下750桁)
4.940656458412465441765687928682213723650598026143247644255856825006755072702087518652998363616359923797965646954457177309266567103559397963987747960107818781263007131903114045278458171678489821036887186360569987307230500063874091535649843873124733972731696151400317153853980741262385655911710266585566867681870395603106249319452715914924553293054565444011274801297099995419319894090804165633245247571478690147267801593552386115501348035264934720193790268107107491703332226844753335720832431936092382893458368060106011506169809753078342277318329247904982524730776375927247874656084778203734469699533647017972677717585125660551199131504891101451037862738167250955837389733598993664809941164205702637090279242767544565229087538682506419718265533447265625e-324
```

# 整数と小数で状況が異なる

小数点以下を含む10進のを F64 に変換したとき、10進数の整数部は2進数に変換しても余程大きな数でない限り同じ整数値ですが、10進小数は2進数表現ができないものが多数あるため一対一対応になりません。

## 整数部の変換

整数部を2進数から10進数への変換の桁数比は概ね $10:3$ で、逆変換も同様です。

| 2進数 | 10進数 |
|--:|--:|
| $11111\ 00111$ | $999$ |
| $11111\ 11111$ | $1,023$ |
| $11110\ 10000\ 10001\ 11111$ | $999,999$ |
| $11111\ 11111\ 11111\ 11111$ | $1,048,575$ |
| $11101\ 11001\ 10101\ 10010\ 01111\ 11111$ | $999,999,999$ |
| $11111\ 11111\ 11111\ 11111\ 11111\ 11111$ | $1,073,741,823$ |
| ... | ... |

```math:
\begin{eqnarray}
\log_{10}{2} & = & 0.30102999566398 \cdots \\
\log_{2}{10} & = & 3.32192809488736 \cdots \\
\end{eqnarray}
```

## 小数部の変換

### 2進数から10進数への変換

小数点部の変換は整数と異なり、2進数から10進数への変換すると桁数比は $1:1$ です。

| 2進数 | 10進数 |
|:--|:--|
| $0.1$ | $0.5$ |
| $0.01$ | $0.25$ |
| $0.001$ | $0.125$ |
| $0.0001$ | $0.0625$ |
| $0.00001$ | $0.03125$ |
| ... | ... |

小数点以下の桁数を $n$ として、小数部の2進数値を $N = 0.\ b_1 b_2 b_3 \cdots b_{n-1}\ 1_n$ から整数 $B$

```math:
B = 2^n \ N = \left( b_1 b_2 b_3 \cdots b_{n-1}\ 1_n \right)_{(2)}
```

とします。式を弄ると

```math:
N = \frac{ B }{ 2^n } = \frac{ B }{ 2^n } \frac{ 5^n }{ 5^n } = \frac{ D = B \times 5^n }{ 10^n }
```

$B$ が2進数表現、$D$ が10進数表現に対応していて、分母の指数から同じ桁数であることが分かります。さらに $D$ は $5^n$ の倍数になります。

### 10進数から2進数への変換

小数点を除いた10進表記を $D$、小数点以下の桁数を $n$ とすると

```math:
N = \frac{ D }{ 10^n }
```

の演算が実行されるでしょう。結果、10進数から2進数への変換では循環小数になるものが多発します。例えば10進数の $0.2$ を2進数にすると

```math:
0.0011\ 0011\ 0011\ 0011 \cdots
```

になります。10進数の $1 \div 3 = 0.3333\cdots$ と同様に $2 \div 10 = 1 \div 5$ の2進数表現では循環小数になります。以下、2進数での計算。

```console:
10進数[ 2] = 2進数[  10]
10進数[10] = 2進数[1010]

         0.0011 0011
      ---------------
 1010 ) 10.000
         1 010
        ------
           1100
           1010
           ----
             10 000
              1 010
             ------
                1100
                1010
                ----
                  10
```

F64 への変換は、手元の環境では

```Python:Python
>>> (2/10).hex()
'0x1.999999999999ap-3'
```

```text:
16進表記: 0x1.999999999999ap-3
厳密値: 0.200000000000000011102230246251565404236316680908203125
```

になります。ゼロ方向丸めの環境では

```text:
16進表記: 0x1.9999999999999p-3
厳密値: 0.1999999999999999833466546306226518936455249786376953125
```

となるかもしれません。$0.2$ は2つの数値の間にあります。

```text:
0x1.9999999999999p-3 = 0.1999999999999999833466546306226518936455249786376953125
0x1.999999999999ap-3 = 0.200000000000000011102230246251565404236316680908203125
```

# F64 の10進数近似

F64 では負符号を省くと、指数部 $E$ (11ビット)と仮数部 $T$ (52ビット)から

```math:
\begin{eqnarray}
n & = & E - 1023 \\
N & = & 2^n \times (0 + \frac{ T }{ 2^{52} }) \ \cdots\ (E = 0) \\
N & = & 2^n \times (1 + \frac{ T }{ 2^{52} }) \ \cdots\ (E \ne 0) \\
\end{eqnarray}
```

の値になります。$T$ が $52$ ビットなので、さらに下位桁の範囲で $N$ には幅(許容誤差)があると考えることができます。次の場合

1. 算出した F64 値を10進数表記
1. 表記された10進数を F64 値として使用

F64 では同値となる範囲で短い10進数表記を探ることができます。しかし、浮動小数点演算プロセッサの丸めや10進数文字列からの F64 値変換処理次第で許容誤差が変わるので、開発環境に合った許容誤差を模索する必要があります。

## サンプルプログラム

```Python:Python
import math
import struct


def altfmtf(value, exact=None, inf='inf', nan='nan'):
    """
    float 型の絶対値10進文字列を得る

    value: 数値
    exact: 厳密値
    inf  : 無限大の文字列
    nan  : 非数の文字列
    """

    LOG_2 = 0.30103  # log10(2) ≒ 0.30103

    B_MAN = 52
    B_EXP = 11

    M_MAN = (1 << B_MAN) - 1
    M_EXP = (1 << B_EXP) - 1

    O_EXP = M_EXP >> 1

    P_MAN = 0
    P_EXP = P_MAN + B_MAN

    B_FX1 = O_EXP + B_MAN - 1

    FX_POW5 = pow(5, B_FX1)
    FX_POW10 = pow(10, B_FX1)

    value = float(value)
    avalue = abs(value)
    fp_data = struct.unpack('<Q', struct.pack('<d', value))[0]
    fp_exp = (fp_data >> P_EXP) & M_EXP
    fp_man = fp_data & M_MAN

    if fp_exp == 0 and fp_man == 0:
        return '0.0'
    if fp_exp == M_EXP:
        return (nan if fp_man else inf)
    if exact and fp_exp > B_FX1:
        return str(int(avalue))

    fp_bin = fp_man | (bool(fp_exp) << 52)
    fx_pos = max(0, fp_exp - 1)

    dn_val = (fp_bin << fx_pos) * FX_POW5
    dn_dlt = (1 << fx_pos) * FX_POW5
    dn_lcn = int((dn_dlt.bit_length() - 1) * LOG_2)
    dn_mod = pow(10, dn_lcn)

    def cl_check(val):  # 再現確認
        try:
            return avalue == val / FX_POW10
        except:
            return False

    # 10進数の下位桁から順に '0' にしつつ再現する値を検索
    cl_fix = [dn_val]
    cl_mod = dn_mod
    while not exact:
        cl_1 = (v - v % cl_mod for v in cl_fix)          # 切り捨て値の表
        cl_2 = sum(([t, t + cl_mod] for t in cl_1), [])  # 切り捨て・切り上げ値の表
        cl_z = list(filter(cl_check, cl_2))              # 再現する値の表
        if not cl_z:
            break
        cl_fix = cl_z
        cl_mod *= 10
    cl_val = sorted(cl_fix, key=lambda v: abs(v - dn_val))[0]  # 元値に最も近い10進数を選択

    s_tmp = str(cl_val)
    s_val = '0' * max(0, 1 + B_FX1 - len(s_tmp)) + s_tmp
    s_int = s_val[:-B_FX1]
    s_dec = s_val[-B_FX1:].rstrip('0')
    return f'{s_int}.{s_dec}'.rstrip('.')


def cvtexp(nstr, inf='inf', nan='nan'):
    """通常10進数を指数形式へ変換"""

    if nstr in (inf, nan):
        return nstr

    s_int, s_dec = (nstr.split('.') + [''])[:2]

    if s_int == '0':
        e_dec = len(s_dec)
        s_dec = s_dec.lstrip('0')
        s_num = f'{s_dec[:1]}.{s_dec[1:]}'.rstrip('.')
        return f'{s_num}e-{e_dec - len(s_dec) + 1:02}'

    e_int = len(s_int)
    s_num = (s_int + s_dec).rstrip('0')
    s_num = f'{s_num[:1]}.{s_num[1:]}'.rstrip('.')
    return f'{s_num}e+{e_int - 1:02}'


def ftopm(value):
    return '-' if value < 0 else ''


def ftoas(value):
    fv = altfmtf(value)
    ev = cvtexp(fv)
    return ftopm(value) + (fv if len(fv) <= len(ev) else ev)

def ftoax(value):
    return ftopm(value) + altfmtf(value, exact=True)


def test(value):
    sv = ftoas(value)
    xv = ftoax(value)
    print(f'16進 : {value.hex()}')
    print(f'厳密 : {xv}')
    print(f'repr : {value}')
    print(f'ftoas: {sv}')
    return
```

```Python:実行結果
Python 3.9.6 (default, Nov 11 2024, 03:15:39)
[Clang 16.0.0 (clang-1600.0.26.6)] on darwin
>>> test(0.2)
16進 : 0x1.999999999999ap-3
厳密 : 0.200000000000000011102230246251565404236316680908203125
repr : 0.2
ftoas: 0.2
>>> test(float.fromhex('0x1.9999999999999p-3'))
16進 : 0x1.9999999999999p-3
厳密 : 0.1999999999999999833466546306226518936455249786376953125
repr : 0.19999999999999998
ftoas: 0.19999999999999998
>>> test(0.1)
16進 : 0x1.999999999999ap-4
厳密 : 0.1000000000000000055511151231257827021181583404541015625
repr : 0.1
ftoas: 0.1
>>> test(float.fromhex('0x1.9999999999999p-4'))
16進 : 0x1.9999999999999p-4
厳密 : 0.09999999999999999167332731531132594682276248931884765625
repr : 0.09999999999999999
ftoas: 0.09999999999999999
>>> test(0.12345678901234567)
16進 : 0x1.f9add3746f65ep-4
厳密 : 0.1234567890123456634920984242853592149913311004638671875
repr : 0.12345678901234566
ftoas: 0.12345678901234566
>>> test(0.123456789012345678)
16進 : 0x1.f9add3746f65fp-4
厳密 : 0.12345678901234567736988623209981597028672695159912109375
repr : 0.12345678901234568
ftoas: 0.12345678901234568
>>> test(1e100)
16進 : 0x1.249ad2594c37dp+332
厳密 : 10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104
repr : 1e+100
ftoas: 1e+100
>>> (1e100).hex()
'0x1.249ad2594c37dp+332'
```

```console:node.jsの場合
Welcome to Node.js v22.13.1
> 0.2
0.2
> 0.1
0.1
> 0.12345678901234567
0.12345678901234566
> 0.123456789012345678
0.12345678901234568
> 1e100
1e+100
```

実行環境は Mac Mini (Intel), macOS 15.2 です。

Python float型の既定の文字列化(repr/strなど)や JavaScript Number 型の既定文字列化も同様の処理が行われているようです。環境違いでも最下位ビット分の誤差が想定できるので、それが許容されるなら、そのまま使ってもよさそうです。(随分前に調べたことを思い出しましたが、すっかり忘れてました...)