浮動小数点数の文字列変換(簡易版)

数学タグが正しいのかどうかわかりませんが…

# 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

# 実装(簡易版)
## 方針
簡易版の方針として、「文字列に変換した値をもう一度`double`に戻した(`strtod`関数とか)とき、元に戻る程度の精度を出す」とします。exactに変換するのもまた楽しいのですが、バッファが無限に(実際のところ3kバイトくらい)膨らみます。
どの程度で切ればよいかというと、$m$の有効ビット数が53であること、したがって有効な数字(digit)が16弱 ( $\log_{10}2^{53}\simeq 15.955$ ) であることから18あればだいたいすみます。
それと、有効bit数が53あるので、これに掛け算割り算を行う都合上64x64=128bitの掛け算と128÷64=64bit(+剰余64bit)という割り算演算が必要です。

## 変数宣言
以下、一時的変数宣言。
```C
    uint64_t mi, mj, h, l;
    int ei;
    unsigned int ex, i, eb, ef;
    size_t pf;
    union {
        double d;
        uint64_t l;
    } d;
    char tmpbuf[20];
```

## ビットの取り出し
まずはビットフィールドを取り出します。`mi`が仮数部、`ei`が指数部になりますが、符号はいきなり返却しても大丈夫。
```C
    d.d = arg;
    mi = d.l & ((1UL << 52) - 1);
    ei = (d.l >> 52) & ((1UL << 11) - 1);
    *sign = d.l >> 63;
```
## 指数仮数の取り出し
例外的なケースを先に対処します。

### e=0
非正規化数($m\neq 0$)の場合、左シフトをかまして$m_{0}=1$になるようにします(下記 $e\neq 0$のケースに合わせてやります)。仮数を左シフトした分指数は引きます($e$が負になります)。

`__builtin_clzl`は左から見て行って0であるビットの数を与えます(0を与えると64を返すと思いきや不定です)。`ei`に対する12という補正は、64ビットに占める仮数でない部分のビット数のことです。
```C
    if (ei == 0) {
        if (mi == 0) {
            *decpt = 0;
            if ((size_t)ndigits > (len - 1))
                ndigits = len - 1;
            memset(buf, '0', ndigits);
            return 0;
        } else {
            eb = __builtin_clzl(mi);
            mi <<= eb - 11;
            ei = - (eb - 12);
        }
    }
```
### e=0x7FF
非数の処置は適当に。バッファの長さより長いコピーはいたしません。
```C
    else if (ei == 0x7ff) {
        *decpt = 0;
        if (len >= 4)
            len = 4;
        if (mi == 0)
            memcpy(buf, "inf", len);
        else
            memcpy(buf, "nan", len);
        return 0;
    }
```
### others
上記以外のケースでは、`mi`に「ケチ表現」で省略されていた最上位ビットを補ってやります($m_0=1$)。
```C
    else {
        /* implicit bit */
        mi |= (1UL << 52);
    }
```
ここまでで`mi`、`ex`に浮動小数点数 `arg`を展開しました。
## 仮数を整数にするような10の乗数を計算
### 指数を求める
$e$の値から、$2^{52} \times 2^e$ が10の何乗くらいなのかを推定します。上記の$m$の設定から、$2^{52}\leq 1m_1m_2...m_{52} < 2^{53}$ですから、もとの数値`arg`が10の何乗程度なのかが推定できます。

以下で、`pf` = 0のときは10の _マイナス_ `ef`乗程度、`pf`が１のときは10の`ef`乗程度となります。`0x4d01`を掛けて16ビット右シフトは、実は$\log_{10}2 = 0.3010$倍です。
```C
    if (ei < 1023 + 52) {
        pf = 0;
        ef = ((((1023 + 52) - ei) * 0x4d01) >> 16) + 1;
    } else {
        pf = 1;
        ef = ((ei - (1023 + 52)) * 0x4d01) >> 16;
    }
```
### 実の値を求める
ここから浮動小数点数計算で10の`ef`乗を求めます。ただし、仮数は`mj`、指数は`ex`に与えられます(53ビットでは精度が足りないので、自力で計算します)。

`tentab`テーブルは論理上は10の$2^{index+1}$乗の値ですが、やや工夫してあって
- 仮数の最上位ビットが1になるように寄せてあるが、小数点の位置は最も左のビットの1つ左
- 指数はこれに合わせてあります

というテーブルです。これに対し
- `mj` の小数点位置は最上位ビットの右

です。64bit×64bitの積をとると、仮数と`mj`の積の小数点の位置は128ビットの最も左側のビットの左に出ますが、これを64ビットの最上位ビットの右に合わせていきます。`+ 1`は切り上げ(の手抜き)です。
```C
    mj = 1UL<<63;
    ex = 0;
    for(i = 0; i < 9; i++) {
        if (ef & (1 << i)) {
            mult128(mj, tentab[i].man, &l, &h);
            if (h & (1UL << 63)) {
                mj = h + 1;
                ex += tentab[i].exp + 1;
            } else {
                mj = ((l >> 63) | (h << 1)) + 1;
                ex += tentab[i].exp;
            }
        }
    }
```
## 引数を整数化
`pf==1`、すなわち`arg`が整数のときは求まった値で`arg`を割ります。`pf==0`のときは`arg`に掛けてやります。これで`mj`が$2^{52}$から$10\times 2^{52}$程度の整数になります。`ex`と`ei`の差が仮数の積を求めるだけではたりない部分で、$2^n$倍する必要があります(下記の補正)。

```C
    if(pf) {
        if (ex != 0) {
/* div */
            mj = div128(mi, 0, mj, NULL);
            pf = ei - ex - 1076;
            mj <<= pf;
        } else {
            mj = mi << (ei - 1075);
        }
        *decpt = ef;
    } else {
/* mult */
        mult128(mi, mj, &l, &h);
        pf = ei + ex - 1074;
        mj = h << pf;
        *decpt = -ef;
    }
```
## 整数を十進変換
`mj`が整数になったので、10進数に変換し、桁数を記録。そのぶん`decpt`を右に移動します。以下のコードは下の桁から、`tmpbuf`の一番後ろから使うように変換していて、出口では`tmpbuf[ei+1]`が0にならないようになっております(バッファの先頭バイトは"0"であってはいけません)。
```C
    for(ei = sizeof(tmpbuf) - 1; mj; ei--) {
        tmpbuf[ei] = '0' + (mj % 10);
        mj /= 10;
    }
    *decpt += sizeof(tmpbuf) - (ei + 1);
```
## あとしまつ
バッファに埋め込みなおして終わりです。
```C
    for(pf = 0, ex = ei + 1; pf < len && pf < (size_t)ndigits; pf++, ex++) {
        if (ex >= sizeof(tmpbuf)) {
           buf[pf] = 0;
        } else {
           buf[pf] = tmpbuf[ex];
        }
    }
    buf[pf] = 0;
    return 0;
```
# 定数表
これを見ると index=5 (10^32) から精度が怪しい感じになっていると思われますが、なんとかセーフのようです。
```C
#include <stdint.h>
static const struct tentab {
    uint64_t man;
    uint16_t exp;
} tentab[] = {
  { 0xa000000000000000UL, 3 },
  { 0xc800000000000000UL, 6 },
  { 0x9c40000000000000UL, 13 },
  { 0xbebc200000000000UL, 26 },
  { 0x8e1bc9bf04000000UL, 53 },
  { 0x9dc5ada82b70b59dUL, 106 },
  { 0xc2781f49ffcfa6d5UL, 212 },
  { 0x93ba47c980e98cdfUL, 425 },
  { 0xaa7eebfb9df9be8dUL, 850 },
  { 0xe319a0aea60e3c87UL, 1700 },
};
```