浮動小数点数の文字列変換(メモリ食い版)

こんどは浮動小数点数をexactに変換していきます。

前半は前の記事( https://qiita.com/nekopon44/items/bfad4e4438787b271d18 )と同じ。
# Cに関するお断り
最近の Cでは共用体(union)を使うときは「代入した型でしか取り出してはいけない」らしいのですが、以下のコードではビットフィールドの取り出しに使っております。

# IEEE754: 準備
binary64形式 (`double`)を扱います。ビット列は$s e_1e_2e_3...e_{11} m_1m_2...m_{52}$と表現できて、
- $s$ -- 符号; 1のときは負の数
- $e = e_1...e_{11}$ -- 指数
- $m = m_1...m_{52}$ -- 仮数

と呼ばれます。$e$ が 0または 0x7FFでないときは、
```math
1.m \times 2^{e - 1022}
```
を表します。ここで出てきた$1.$は「隠れた1」「ケチ表現」などと呼ばれるもので、有効ビット数を稼ぐために生まれた何かであります。また1022は指数部を負にしないためのバイアスと呼ばれる数値です。
$e=0$のときは
```math
0.m \times 2^{-1022}
```
を表します。ケチ表現がなくなっているのと、バイアスが「1023になってない」ところが大事なポイントです。$e=0、m=0$であればふつうに0です。
$e=\mathrm{0x7FF}$のときは数値とは取り扱われず、
- $m=0 : \infty$
- $m\neq 0 :$ Not a number (NaN)

と表現されます。

`long double`形式は実はポータブルでなく、分解はコンパイラ依存になってきます。まあ展開ができてビット数がそれなりな(64bitではちょっと足りません)乗除算ができれば、以下と同じ仕掛けで`qecvt_r`を書くことはできるはずであります。
# ecvt_r
Deprecatedな関数ですが、これを実装していきます。現代的な手法は`sprintf`を使うこと。
```C
int ecvt_r(double arg, int ndigits, int *decpt, int *sign, char *buf, size_t len);
```
引数
- `arg` : 変換する値
- `ndigits` : 変換する桁数
- `decpt` : [返却] 小数点の位置。0のときバッファの先頭
- `sign` : [返却] 負のときは0でない値
- `buf`, `len` : バッファ

返却値
- 正常に変換できれば0

# 実装(メモリ食い版)
## 方針
メモリ食い版は、以下の仕掛けを使います。
- まず整数部と小数部に分解
- これを多倍長演算を使って9桁ごとに変換

多倍長演算は、複数の変数をつかって1つの数値を表現する仕組みです。意味としては進数表示と同じことをやるのですが、配列`m`が表す数値を$\sum_{k=0}^{n}m_k\times B^k$、ないし$\sum_{k=0}^{n}m_k\times B^{-k}$とします。ここで定数$B$は「基数」と呼ばれる数で、10進数であれば10になるような数です。
演算の数をなるべく少な目にすることを考えて、配列の1エントリの大きさを32bit`uint32_t`としたので基数は$2^{32}$です。

そして`uint32_t`にはいる10の累乗の数値で一番大きいのが$10^9$です($10^9 < 2^{32}$)。ということで以下の`__cvt_int`、`__cvt_frac`の多倍長バッファでは、一度に10進数9桁を変換していきます。

## 整数部の変換
```C
static int __cvt_int(uint64_t man, int exp, uint32_t * intbuf);
```
引数
- `man`: 整数の仮数部
- `exp`: 整数の指数部
- `intbuf`: 返り値のバッファ

返却値は使った`intbuf`バッファの数です。引数たる整数は $man\times 2^{exp}$ です。
### 実装
1024ビットが入るような配列を用意し、これを多倍長演算に使います。
```C
    uint32_t a[(DBL_MAX_EXP + DBL_MANT_DIG + 31)/32];   /** Real bit array */
    uint64_t ml;
    int i, j, c, k;

    /* clear temporal storage */
    memset(a, 0, sizeof(a));
```
`man`と`exp`を配列`a`に転写していきます。コメントにあるように、indexが0のときがもっとも下位です。`j`は配列`a`で使っているところの最大値、言い換えると`a[j]`までは有効な値が入っている、ということです。

後半みるとわかるとおり、最大3つのエントリにまで`man`の値が入ります。
```C
#define DBL_MANT_DIG_1  (DBL_MANT_DIG-1)
    /**
     * load "bit array". Least significant word is a[0]
     */
    if (exp < DBL_MANT_DIG_1) {
        if ((a[1] = (uint32_t) (man >> (64 - DBL_MANT_DIG_1 - exp))) != 0)
            j = 1;
        else
            j = 0;
        a[0] = (uint32_t) (man >> (DBL_MANT_DIG_1 - exp));
    } else {
        j = (exp - DBL_MANT_DIG_1) / 32;
        k = (exp - DBL_MANT_DIG_1) % 32;
        if (k == 0) {
            a[j + 1] = (uint32_t) (man >> 32);
            a[j] = (uint32_t) man;
            j += 1;
        } else {
            a[j + 2] = (uint32_t) (man >> (64 - k));
            a[j + 1] = (uint32_t) (man >> (32 - k));
            a[j] = (uint32_t) (man << k);
            j += 2;
        }
    }
```

$10^9$で割った商とあまりを求めます。剰余は10進9桁の変換に相当します。…というのを、商が0になるまで繰り返してやります。
```C
    /**
     * divmod by 1e9
     */
    c = 0;
    do {
        ml = 0;
```
割り算は上の桁から。(`j`は最上位のエントリを指しています)
```C
        for (i = j; i >= 0; i--) {
```
上位のエントリの剰余と現在のエントリを足してやります。基数$2^{32}$ですから、上位のエントリ分は左32ビットシフト。また計算は64ビットが必要です。
得られた値を$10^9$で割って商と剰余を求めます。商は次回のループで使い、剰余は下位の桁へ伝達します。
```C
            ml = (ml << 32) + a[i];
            a[i] = (uint32_t) (ml / 1000000000UL);
            ml = ml - a[i] * 1000000000UL;
```
最下位のエントリまで計算したら、剰余(`ml`)が9桁分の変換値になりますので、`intbuf`へ1桁記録。
```C
        }
        intbuf[c++] = (uint32_t) ml;
```
最上位のエントリが0になっていたら、そのエントリはもう計算する必要はありません。そして全部のエントリが吐き出されたら計算自体は終了です。
```C
        if (a[j] == 0) {
            --j;
        }
    } while (j >= 0);
```
最後にちょっとした補正として、`intbuf`の最上位が0になる場合は、そのエントリを削ってしまいます。
```C

    /**
     * To workaround the most significant word may be zero
     */
    if (intbuf[c - 1] == 0)
        c--;
    return c;
}
```
## 小数部の変換
上とほぼ同じ技術です。

配列`a`の大きさは1074ビット(1022+52ビット)必要となります。
```C
static int __cvt_frac(uint64_t man, int exp, uint32_t * buf)
{
    uint32_t a[(-DBL_MIN_EXP + DBL_MANT_DIG + 31) / 32];
    int j, k, fl, t;
    uint64_t ml, mk;

    memset(a, 0, sizeof(a));
```

さっきと同じような仕掛けで、配列`a`を埋めていきます。このケースでは、数値は$\sum_{i=0}^{j}a_i\times 2^{-32\times(i+1)}$となっていて、小数点の位置は`a[0]`の最上位ビットの左です。
```C
    /**
     * load "bit array". Most significant word is a[0]
     */
    if (exp < 64 - DBL_MANT_DIG_1) {
        a[0] = (uint32_t) (man >> (exp + (DBL_MANT_DIG_1 - 32)));
        if ((a[1] = (uint32_t) (man << (DBL_MANT_DIG_1 - exp))) != 0)
            j = 1;
        else
            j = 0;
    } else {
        exp -= 64 - DBL_MANT_DIG_1;
        k = exp % 32;
        j = exp / 32;
        if (k == 0) {
            a[j] = (uint32_t) (man >> 32);
            if ((a[j + 1] = (uint32_t) (man >> 0)) != 0)
                j++;
        } else { 
            fl = 0;
            a[j] = (uint32_t) (man >> (32 + k));
            if ((a[j + 1] = (uint32_t) (man >> k)) != 0)
                fl = 1;
            if ((a[j + 2] = (uint32_t) (man << (32 - k))) != 0)
                fl = 2;
            j += fl;
        }
    }
```
こんどは$10^9$を掛けていきます。小数部に$10^9$を掛けてやると、$10^9-1$までの整数が得られるはずです。
```C
    /** Generate digits. Mult 1e9 so we can 9 digits at once */
    t = 0;
    while (j > 0 || a[0] != 0) {
        ml = 0;
```
掛け算は下位のエントリからです。
```C
        for (k = j; k >= 0; k--) {
```
32bit×32bit=64bitの積を得て、上位32ビットは上のエントリへもちこしていきます。
```C
            mk = a[k];
            mk = mk * 1000000000UL;
            ml = (ml >> 32) + mk;
            a[k] = (uint32_t) ml;
```
`a[0]`の$10^9$倍した持ち越しぶんこそ変換結果です。
```C
        }
        buf[t++] = (uint32_t) (ml >> 32);
```
最下位のエントリが0になっていたら、そのエントリは計算する必要がなくなります。
```C
        if (a[j] == 0)
            j--;
```
変換したエントリ数を返して終わりです。
```C
    }
    return t;
}
```
## 分解ととりまとめ
簡易版と同じように前処理しておきます(細かい説明はパス)。`ibuf`の大きさはやや大きめに見積もっております(1075桁/9桁は必要)。
```C
static const char tostr[] = "0123456789";

int
ecvt_r(double arg, int ndigits, int *decpt, int *sign, char *buf, size_t len)
{
    uint64_t man, mani, manf;
    int exp, expf;
    union {
        uint64_t l;
        double d;
    } d;
    int i, j, ilen;
    unsigned int fl, nd;
    size_t icnt;
    uint32_t ibuf[128], tbuf;

    if (!decpt || !sign || !buf || len < 2 || ndigits <= 0
        || (unsigned) ndigits >= len)
        return -1;

    nd = ndigits;
    d.d = arg;
    *sign = d.l >> 63;
    man = d.l & ((1UL << 52) - 1);
    exp = (d.l >> 52) & 0x7ff;
    if (exp == 0 && man == 0) {
        /* fully zero */
        memset(buf, tostr[0], nd);
        *decpt = 1;
        return 0;
    }
    if (exp == 0x7ff) {
        /* inf/nan */
        if (man) {
            strncpy(buf, "nan", len);
        } else {
            strncpy(buf, "inf", len);
        }
        return 0;
    }

    /* nonzero, may denormalized */
    if (exp != 0)
        man |= (1UL << 52);
```
ここまでで`man`に仮数、`exp`に指数+バイアスが得られます。

バイアス付き指数が1023以上(指数が0以上)のときは、整数部があります。たとえば、$e=1023$のときは、ちょうど「隠れた1」が整数になります。小数部はマスクします。
```C
    if (exp >= 1023) {          /* has integer part */
        if (exp < 1023 + 52) {
            mani = man & ~((1UL << (1023 + 52 - exp)) - 1);
        } else {
            mani = man;
        }
    } else {
        mani = 0;
    }
```

 以下では小数部を取り出しています。$e$が1023以上の時は整数部がありますから、先ほどとは逆に左側のビットをマスクします。
```C
    if (exp < 1023 + 52) {      /* has fractal part */
        if (exp >= 1023) {      /* has integer part */
            manf = man & ((1UL << (1023 + 52 - exp)) - 1);
```
以下の計算は`manf`の最上位ビットをなるべく左に寄せています(その分指数を調整する)。
```C
            expf = __builtin_clzl(1) - __builtin_clzl(manf);
            manf <<= 52 - expf;
            expf = 1023 - exp + 52 - expf;
        } else if (exp != 0) {  /* normalized */
            manf = man;
            expf = 1023 - exp;
        } else {                /* denormalized */
            manf = man;
            expf = 1022;
        }
    } else {
        manf = 0;
    }
```
整数部を変換します。
```C
    if (mani != 0) {
        ilen = __cvt_int(mani, exp - 1023, ibuf);
```
ひどいコードですが、最上位エントリの桁数を調べます。最上位エントリのぶんがわかれば、残りはすべて9桁です。
```C
        ilen--;
        if (ibuf[ilen] >= 100000000U)
            fl = 9;
        else if (ibuf[ilen] >= 10000000U)
            fl = 8;
        else if (ibuf[ilen] >= 1000000U)
            fl = 7;
        else if (ibuf[ilen] >= 100000U)
            fl = 6;
        else if (ibuf[ilen] >= 10000U)
            fl = 5;
        else if (ibuf[ilen] >= 1000U)
            fl = 4;
        else if (ibuf[ilen] >= 100U)
            fl = 3;
        else if (ibuf[ilen] >= 10U)
            fl = 2;
        else
            fl = 1;
        *decpt = ilen * 9 + fl;
```
今度は10進1桁ごとの形に変換します。

ここで、結果のバッファ`buf`がたりるかどうかを`ndigits`の値でチェックしていきます。バッファはいっぱいまで使います。
```C
        for (icnt = 0; ilen >= 0 && icnt <= nd; ilen--, fl = 9) {
            tbuf = ibuf[ilen];
            for (i = fl - 1; i >= 0; i--) {
                if (icnt + i <= nd) {
                    buf[icnt + i] = tostr[tbuf % 10];
                }
                tbuf /= 10;
            }
            icnt += fl;
        }
```
`icnt`は整数の桁数となります。もし最終桁が書かれていなかったら"0"にしておきます。
```C
        /* fixup length */
        if (icnt > nd + 1) {
            icnt = nd + 1;
        }
        if (nd > icnt) {
            buf[icnt] = tostr[0];
        }
```
`mani=0`のケースは整数部が無い場合ですね。
```C
    } else {
        *decpt = 0;
        icnt = 0;
        buf[0] = tostr[0];      /* may "round up handler" use it */
    }
```
まだ桁が余っているかどうかを確認。バッファに余裕があり、かつ小数部が0でないのなら、小数を変換していきます。やりかたは先ほどと同じ。
```C
    if (nd >= icnt) {
        if (manf != 0) {
            ilen = __cvt_frac(manf, expf, ibuf);
```
ただし、「整数部が0の場合に限り」小数部の最上位エントリの桁数が必要です。バッファの最上位バイトが"0"であってはいけないので、このぶんを変換せずに取り除いてやらねばならないのです。
```C
            if (icnt == 0) {
                /* eliminate leading zeros */
                for (i = 0; ibuf[i] == 0; i++);

                if (ibuf[i] >= 100000000U)
                    fl = 9;
                else if (ibuf[i] >= 10000000U)
                    fl = 8;
                else if (ibuf[i] >= 1000000U)
                    fl = 7;
                else if (ibuf[i] >= 100000U)
                    fl = 6;
                else if (ibuf[i] >= 10000U)
                    fl = 5;
                else if (ibuf[i] >= 1000U)
                    fl = 4;
                else if (ibuf[i] >= 100U)
                    fl = 3;
                else if (ibuf[i] >= 10U)
                    fl = 2;
                else
                    fl = 1;
                *decpt = fl - ((i + 1) * 9);
            } else {
                /* dont eliminate leading zeros */
                i = 0;
                fl = 9;
            }
```
`fl`が変換すべき桁数となっていますので、エントリを10進数に変換していきます。

```C
            for (; i < ilen && icnt <= nd; i++, fl = 9) {
                tbuf = ibuf[i];
                for (j = fl - 1; j >= 0; j--) {
                    if (icnt + j <= nd) {
                        buf[icnt + j] = tostr[tbuf % 10];
                    }
                    tbuf /= 10;
                }
                if (nd - icnt >= fl) {
                    icnt += fl;
                } else {
                    icnt += nd - icnt + 1;
                }
            }
        } else {
            while (icnt <= nd + 1) {
                buf[icnt++] = tostr[0];
            }
        }
    }
 
```
上記で見たとおり、必要な桁数(`ndigit`)より1桁余分に変換しております。これを使って、最下位の四捨五入を実行します。
```C
    icnt--;
    /* round up */
    if (buf[icnt] >= tostr[5]) {
        for (i = icnt - 1; i >= 0; i--) {
            if (++buf[i] <= tostr[9]) {
                break;
            }
            buf[i] = tostr[0];
        }
        /* overflow */
        if (i < 0) {
            buf[0] = tostr[1];
            (*decpt)++;
        }
    }
```
ラップアラウンド処理の後をもみ消して、okを返してやります。
```C
    buf[icnt] = 0;
    return 0;
}
```