こんどは浮動小数点数を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でないときは、
1.m \times 2^{e - 1022}
を表します。ここで出てきた$1.$は「隠れた1」「ケチ表現」などと呼ばれるもので、有効ビット数を稼ぐために生まれた何かであります。また1022は指数部を負にしないためのバイアスと呼ばれる数値です。
$e=0$のときは
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を使うこと。
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エントリの大きさを32bituint32_tとしたので基数は$2^{32}$です。
そしてuint32_tにはいる10の累乗の数値で一番大きいのが$10^9$です($10^9 < 2^{32}$)。ということで以下の__cvt_int、__cvt_fracの多倍長バッファでは、一度に10進数9桁を変換していきます。
整数部の変換
static int __cvt_int(uint64_t man, int exp, uint32_t * intbuf);
引数
-
man: 整数の仮数部 -
exp: 整数の指数部 -
intbuf: 返り値のバッファ
返却値は使ったintbufバッファの数です。引数たる整数は $man\times 2^{exp}$ です。
実装
1024ビットが入るような配列を用意し、これを多倍長演算に使います。
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の値が入ります。
#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になるまで繰り返してやります。
/**
* divmod by 1e9
*/
c = 0;
do {
ml = 0;
割り算は上の桁から。(jは最上位のエントリを指しています)
for (i = j; i >= 0; i--) {
上位のエントリの剰余と現在のエントリを足してやります。基数$2^{32}$ですから、上位のエントリ分は左32ビットシフト。また計算は64ビットが必要です。
得られた値を$10^9$で割って商と剰余を求めます。商は次回のループで使い、剰余は下位の桁へ伝達します。
ml = (ml << 32) + a[i];
a[i] = (uint32_t) (ml / 1000000000UL);
ml = ml - a[i] * 1000000000UL;
最下位のエントリまで計算したら、剰余(ml)が9桁分の変換値になりますので、intbufへ1桁記録。
}
intbuf[c++] = (uint32_t) ml;
最上位のエントリが0になっていたら、そのエントリはもう計算する必要はありません。そして全部のエントリが吐き出されたら計算自体は終了です。
if (a[j] == 0) {
--j;
}
} while (j >= 0);
最後にちょっとした補正として、intbufの最上位が0になる場合は、そのエントリを削ってしまいます。
/**
* To workaround the most significant word may be zero
*/
if (intbuf[c - 1] == 0)
c--;
return c;
}
小数部の変換
上とほぼ同じ技術です。
配列aの大きさは1074ビット(1022+52ビット)必要となります。
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]の最上位ビットの左です。
/**
* 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$までの整数が得られるはずです。
/** Generate digits. Mult 1e9 so we can 9 digits at once */
t = 0;
while (j > 0 || a[0] != 0) {
ml = 0;
掛け算は下位のエントリからです。
for (k = j; k >= 0; k--) {
32bit×32bit=64bitの積を得て、上位32ビットは上のエントリへもちこしていきます。
mk = a[k];
mk = mk * 1000000000UL;
ml = (ml >> 32) + mk;
a[k] = (uint32_t) ml;
a[0]の$10^9$倍した持ち越しぶんこそ変換結果です。
}
buf[t++] = (uint32_t) (ml >> 32);
最下位のエントリが0になっていたら、そのエントリは計算する必要がなくなります。
if (a[j] == 0)
j--;
変換したエントリ数を返して終わりです。
}
return t;
}
分解ととりまとめ
簡易版と同じように前処理しておきます(細かい説明はパス)。ibufの大きさはやや大きめに見積もっております(1075桁/9桁は必要)。
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」が整数になります。小数部はマスクします。
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以上の時は整数部がありますから、先ほどとは逆に左側のビットをマスクします。
if (exp < 1023 + 52) { /* has fractal part */
if (exp >= 1023) { /* has integer part */
manf = man & ((1UL << (1023 + 52 - exp)) - 1);
以下の計算はmanfの最上位ビットをなるべく左に寄せています(その分指数を調整する)。
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;
}
整数部を変換します。
if (mani != 0) {
ilen = __cvt_int(mani, exp - 1023, ibuf);
ひどいコードですが、最上位エントリの桁数を調べます。最上位エントリのぶんがわかれば、残りはすべて9桁です。
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の値でチェックしていきます。バッファはいっぱいまで使います。
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"にしておきます。
/* fixup length */
if (icnt > nd + 1) {
icnt = nd + 1;
}
if (nd > icnt) {
buf[icnt] = tostr[0];
}
mani=0のケースは整数部が無い場合ですね。
} else {
*decpt = 0;
icnt = 0;
buf[0] = tostr[0]; /* may "round up handler" use it */
}
まだ桁が余っているかどうかを確認。バッファに余裕があり、かつ小数部が0でないのなら、小数を変換していきます。やりかたは先ほどと同じ。
if (nd >= icnt) {
if (manf != 0) {
ilen = __cvt_frac(manf, expf, ibuf);
ただし、「整数部が0の場合に限り」小数部の最上位エントリの桁数が必要です。バッファの最上位バイトが"0"であってはいけないので、このぶんを変換せずに取り除いてやらねばならないのです。
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進数に変換していきます。
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桁余分に変換しております。これを使って、最下位の四捨五入を実行します。
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を返してやります。
buf[icnt] = 0;
return 0;
}