数学タグが正しいのかどうかわかりませんが…
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
実装(簡易版)
方針
簡易版の方針として、「文字列に変換した値をもう一度doubleに戻した(strtod関数とか)とき、元に戻る程度の精度を出す」とします。exactに変換するのもまた楽しいのですが、バッファが無限に(実際のところ3kバイトくらい)膨らみます。
どの程度で切ればよいかというと、$m$の有効ビット数が53であること、したがって有効な数字(digit)が16弱 ( $\log_{10}2^{53}\simeq 15.955$ ) であることから18あればだいたいすみます。
それと、有効bit数が53あるので、これに掛け算割り算を行う都合上64x64=128bitの掛け算と128÷64=64bit(+剰余64bit)という割り算演算が必要です。
変数宣言
以下、一時的変数宣言。
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が指数部になりますが、符号はいきなり返却しても大丈夫。
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ビットに占める仮数でない部分のビット数のことです。
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
非数の処置は適当に。バッファの長さより長いコピーはいたしません。
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$)。
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$倍です。
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は切り上げ(の手抜き)です。
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$倍する必要があります(下記の補正)。
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"であってはいけません)。
for(ei = sizeof(tmpbuf) - 1; mj; ei--) {
tmpbuf[ei] = '0' + (mj % 10);
mj /= 10;
}
*decpt += sizeof(tmpbuf) - (ei + 1);
あとしまつ
バッファに埋め込みなおして終わりです。
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) から精度が怪しい感じになっていると思われますが、なんとかセーフのようです。
#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 },
};