What goes on/what's going on:2021年04月16日分

@hdtoaでprintf(3)の"%a"書式指定子をエミュってみる。

これまでe/f/g変換は10進表記だったのだけれど、今回できるかなのa変換は16進表記でありP表記(P notation)に変換する。

んでP表記(P-notation)ってのはこんなん。

(?<符号>[+\-])?(?<仮数部>0[xX](?<整数部>[[:xdigit:]])((?<小数点>\.)(?<小数部>[[:xdigit:]]+))?)(?<p表記>[pP])(?<指数部>(?<指数符号>[+\-])(?<指数>[[:digit:]]+))

はい、E表記と比べると仮数部が0xh.hhhhで指数部はp±d形式になっとります。 そして指数はE表記では2桁以上だったのが1桁以上でゼロ埋め不要になっとります。

そんでP表記でprintf(3)なら"%p"だろ！となりそうなもんだけど、残念ながらとっくの昔に"%p"書式指定子はポインタを表示するものとして先約済なので、かわりに"%a"書式指定子を使うとC99で定められたのですわ。 なぜ数ある中から"%a"が選ばれたのか残念ながらワイ経緯知らんしそれっぽい理由も思いつかないし調べる気力も無いので以下略

前回例にとりあげた0.1(これ自体循環小数だけど)から1bit前後にずらした値は10進だと精度によってはどれも0.1に丸められてしまったけれど

$ cat >unko.c
#include <stdio.h>
#include <string.h>

#define arraycount(array)       (sizeof(array)/sizeof(array[0]))

union double2bytes {
	double dvalue;
	char bvalue[8];
};

int
main(void)
{
	const unsigned char *s[] = {
		"\x99\x99\x99\x99\x99\x99\xb9\x3f",
		"\x9a\x99\x99\x99\x99\x99\xb9\x3f",
		"\x9b\x99\x99\x99\x99\x99\xb9\x3f",
	};
	union double2bytes u;
	size_t i;

	for (i = 0; i < arraycount(s); ++i) {
		memcpy(u.bvalue, s[i], sizeof(u.bvalue));
		printf("dvalue[%g]\n", u.dvalue);
	}
}
^D
$ make unko
$ ./unko
dvalue[0x1.9999999999999p-4]
dvalue[0x1.999999999999ap-4]
dvalue[0x1.999999999999bp-4]

はい、16は2のべき乗なので正確に表現できとりますね、そこの36bitワードマシンとか二進化十進数に一言お持ちのお方、ご静粛にご着席ください。

ちなみにこの仮数部の1bitの差(つまり1と1より仮数部1bitだけ大きい数字の差)をマシンイプシロン(青い装甲騎兵ではない)あるいは計算機イプシロンと呼ぶ(doubleだとほとんどの実装において「0x1p-52」となる)んだけどその差を誤差無く表現できるってわけよ。

どうでもいいけどマシンイプシロンについてはCの規格でfloat.hにDBL_EPSILONという定数を用意しろとある、Nだとfloat.hからインクルードされてるsys/float_ieee754.hで

#if __STDC_VERSION__ >= 199901L
#define DBL_EPSILON     0x1.0p-52
...
#else
#define DBL_EPSILON     2.2204460492503131E-16
...
#endif

となっとりますな、なお環境によってはコンパイラの事前定義マクロ__DBL_EPSILON__を使うのでそっちは以下で確認して。

$ echo __DBL_EPSILON__ | gcc -E -
# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "<stdin>"
((double)2.22044604925031308084726333618164062e-16L)

@まさかの時の現代貨幣理論

経済学部で音楽専攻してたワイはギターは絶対にインチネジ派なので度量法をメートルで統一とか絶対に許さないよ…なのだが、貨幣はもう16進にして

• 0x186a0p+0円記念硬貨
• 0x2710p+0円札
• 0x1388p+0円札
• 0x7d0p+0円札
• 0x3e8p+0円札
• 0x1f4p+0円硬貨
• 0x64p+0円硬貨
• 0x32p+0円硬貨
• 0xap+0円硬貨

を発行すればなにかと開発力の低さが社会信用問題になるジャパニーズ金融機関もシステムトラブル発生しなくなるというご提案、いかがですかFinTechとかいう横文字ギョーカイや財務省官僚の方々。

これが実現すれば、右や左の旦那さま恵まれぬ私にどうか0x1.388p+13円お恵みを～と物乞ひすれば、1.4円くらいの小銭と勘違いして銀行口座に気軽に振り込んでくれるだろうしwin-winでなのである。

@hdtoaに渡すパラメーター

話を元に戻して"%a"の実装なんだけど、doubleの内部表現を共用体使って引き摺り出せるならP表記への変換なんて自前で簡単に実装できる気がしてくるけど、精度指定した場合有効桁数未満での丸めはやっぱり必要になるのでgdtoaにお任せした方がいい。

そこでご用意されているのがdtoaの16進バージョンであるhdtoa、今回もこいつに頼ることにする。

char *
hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
	    char **rve)
{

dtoaと比べてだいたい引数は一緒なんだけど、差異としては

• dtoaのような変換モード指定はない(そもそも要らない)
• dtoaと違って変換結果を大文字と小文字のどちらにするかの指定は必要になるけど、第二引数に以下のいづれかをを渡すことで実現
  • "0123456789ABCDEF"
  • "0123456789abcdef"

ちゅーあたりか。

ここで一点注意、さっきの俺の口座に0x1.388p+13円を振り込むには0x2710p+0円札プラス手数料が必要と気づいた人はもうお判りだろうけど、E表記と同じくP表記は同じ値であっても複数の表現を持つのだ。

0x1.3p3
0x2.6p2
0x4.cp1

これ仕様では

there is one hexadecimal digit (which shall be non-zero if the argument is a normalized floating-point number and is otherwise unspecified) before the decimal-point character

と小数点の前は一桁だよだけしか決まってないので、どの表現になるか決まらないのだ(%eの場合はこれだけで決定するんだが)。よってどの表現がでてきても正しい実装依存ガチャなのだ。

このことはhdtoaのソースにも書いてあって、例もそこから抜粋したのだが

 * Note that the C99 standard does not specify what the leading digit
 * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
 * as 0x2.6p2 is the same as 0x4.cp3.  This implementation chooses the
 * first digit so that subsequent digits are aligned on nibble
 * boundaries (before rounding).

はい0x4.cp「1」を0x4.cp「3」とtypoしとりますな、やはり県北の土手でプルリク盛りあおうや…

よって%aを実装しても正しく変換できたかのテストケースを書こうしても移植性考えるとちょっとめんどくさい、結果をstrtod(3)で再度doubleに戻した後はコンパイラが対応してりゃP表記リテラルで書いた期待値と比較すればいいけど、未対応だと期待値をバイナリで書かないとならん。まぁ今時さすがに無いとは思いたい。

@hdtoaが返す文字列(head)のフォーマット、そしてexpで返される値とは

とりあえず適当な値入れて呼んでみる。

$ cat >unko.c
#include <stdio.h>
#include <stdlib.h>
#include "gdtoa.h"

int
main(void)
{
	static const char *xdigit = "0123456789abcdef";
	char *head, *tail;
	int exp, neg;

	head = hdtoa(0x1p-52, xdigit, 0, &exp, &neg, &tail);
	if (head == NULL)
		abort();
	printf("head:[%s], exp:[%d], neg:[%d]\n", head, exp, neg);
}
^D
$ make unko
$ ./unko
head:[1], exp:[-51], neg:[0]

はい指数が-51で帰ってきてるのでhdtoaが変換した結果は「0x0.1p-51」になるということやね。

では精度を指定すると？

$ cat >unko.c
#include <stdio.h>
#include <stdlib.h>
#include "gdtoa.h"

int
main(void)
{
	static const char *xdigit = "0123456789abcdef";
	char *head, *tail;
	int exp, neg;

	head = hdtoa(0x1p-52, xdigit, 6, &exp, &neg, &tail);
	if (head == NULL)
		abort();
	printf("head:[%s], exp:[%d], neg:[%d]\n", head, exp, neg);
}
^D
$ make unko
$ ./unko
head:[100000], exp:[-51], neg:[0]

ほうdtoaの時とは違ってhdtoaは右側のゼロ出力するやんけ、「0x0.100000p-51」ってことやな。

(追記) ただし0.0を変換した場合は「000000」にはならず「0」になるのでその時だけは自分でゼロ埋め必要、めんどくさ。

ということでCの仕様である

• 小数点の前に1桁必要
• デフォルトの精度は6ではなく変換された全ての桁

よって「0x1p-52」あるい「0x1.000000p-52」と表示させるにはprecとexpの値の調整が必要になるわけだ。e変換とおんなじだな。

@いつものクッソ適当な実装例

つーことで前々回のコードをちょいちょい弄ってこうなる。

#include <assert.h>
#include <limits.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>

#define PRECISION       0
#define MAXEXPSIZ       3

static const char *xdigit = "0123456789abcdef";

int
main(int argc, char *argv[])
{
	double dvalue = ...
	char *head, *tail, expstr[MAXEXPSIZ + 2];
	int prec = PRECISION, exp, neg, len, zeros, expsig;
	FILE *fp = stdout;

	if (prec >= 0)
		++prec;
	head = hdtoa(dvalue, xdigit, prec, &exp, &neg, &tail);
	len = tail - head;
	if (neg && putc('-', fp) == EOF)
		abort();
	if (exp == INT_MAX) {
		if (fwrite(head, 1, len, fp) != len)
			abort();
	} else {
		if (fwrite("0x", 1, 2, fp) != 2)
			abort();
		if (putc(*head, fp) == EOF)
			abort();
		if (--len > 0) {
			if (putc('.', fp) == EOF)
				abort();
			if (fwrite(&head[1], 1, len, fp) != len)
				abort();
		}
		if (--exp >= 0) {
			assert(exp <= DBL_MAX_10_EXP);
			expsig = '+';
                } else {
			assert(exp >= DBL_MIN_10_EXP);
			expsig = '-';
			exp = -exp;
		}
		tail = &expstr[sizeof(expstr)];
		while (exp > 9) {
			*--tail = '0' + (exp % 10);
			exp /= 10;
		}
		*--tail = '0' + exp;
		*--tail = expsig;
		*--tail = 'p';
		len = &expstr[sizeof(expstr)] - tail;
		if (fwrite(tail, 1, len, fp) != len)
			abort();
	}
	freedtoa(head);
	exit(EXIT_SUCCESS);
}

もう勘所はわかってるのでちょっと考えれば書けますやね。

一点だけ注意、MAXEXPSIZという謎の定数はe変換の時は4でその理由として

DBL_MIN_EXP	-1021
DBL_MAX_EXP	+1024

を格納できうるサイズという説明をした。

今回a変換の場合は3になっとるのだけど、これは同様にC99の5.2.4.1にてIEEE 754を満たしてるなら指数の上限加減が-307～308が最小値となるとあるからの数字。

DBL_MIN_10_EXP	-307
DBL_MAX_10_EXP	308

まぁ念のため実際の値を自分の使ってる環境値を確認しといてください。

(追記) INF/NANの場合dtoaだと9999がexpに返るけど、hdtoaの場合はINT_MAXになるのを忘れててそのままの実装をお出ししてしまったので修正したよ。

@次回予告

とりあえずこれで やる気のない駆け足だけれどもNなんかのprintf(3)のコードを読む際にわけわかんねーとなりがちな浮動小数点数とgdtoaの関係についての解説はおしまい。 なーんだprintf(3)の実装なんて簡単そうと思っていただけたら幸いであるが、お前の書く文章は読みづらいということであればここは俺のチラシの裏なので苦情は/dev/nullへ。

なお次回からはLC_NUMERICによる国際化という地獄につきあってもらう、あと機会があれば"%1$f"みたいなpositional orderもだな…