Not only is the Internet dead, it's starting to smell really bad.:2021年04月15日分

2021/04/15(Thu)

[プログラミング] gdtoaを使う(その4)

前回の続き、e変換のコードも整理して ここに貼っておいた。

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

今回の%gというのは

  • %eとして変換を行い、指数が-4より大きくかつ精度で指定した有効桁数以下なら
    • %fと同じddd.ddd形式で表示する
    • それ以外なら%eと同じd.ddde±dd形式で
  • 精度は%eとは異なり小数点以下の有効桁数ではなく整数部の1桁目も含む
    • なので0を指定した場合は1と解釈される
  • 精度はあくまで%eとして変換するときのもので%f形式の有効桁数は以下にも以上にもなる、よって
    • 有効桁数未満でも(#を指定しない限り)ゼロ埋めしない
    • 〃より大きくとも切り詰めたり丸めたりせず、全ての桁を表示する

という変換を行う、なんかややこしそうに感じるけど前々回と前回の内容を理解してればささっと書ける。

ちなみに一時配布元のねっとり部 *1ことnetlib.orgのページにはgdtoa作者によるg_fmt.cというソースコードがあるのでこの通りに書けばいい勘違いしてしまうのだが、実はこいつはちょっとprintf(3)の%gとは挙動が違う。

 char *
g_fmt(register char *b, double x)
{
	...
	s = s0 = dtoa(x, 0, 0, &decpt, &sign, &se);

はい、モードと有効桁数どちらもゼロ指定で変換しとりますやね。

モード0ってのはこれ。

                0 ==> shortest string that yields d when read in
                        and rounded to nearest.

超訳すると

	0 ==> doubleを読み込んで誤差を最近接丸めした値を最短の文字列にする。

ですかね。

最短の文字ってどゆこと?なんだけど、2進数で実装されてるdoubleで10進数を表現しようとすると割り切れず循環小数になる場合がある、例えば0.1は正確に0.1ではなくクッソ適当に精度30とか指定すると

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

int
main(int argc, char *argv[])
{
	printf("%.30f\n", 0.1);
}
^D
$ make unko
$ ./unko
0.100000000000000005551115123126

ちゅーかんじで誤差が存在ってのが判ると思う。

しかしdtoaにモード0を指定すると有効桁数に30とか指定しても

$ cat >unko.c
#include <stdio.h>
#include "gdtoa.h"
int
main(void)
{
	char *head, *tail;
	int exp, neg;

	head = dtoa(0.1, 0, 30, &exp, &neg, &tail);
	printf("head:[%s], exp:[%d\]n", head, exp);
}
^D
$ make unko
$ ./unko
head:[1], exp:[0]

このように無視される、どうやら実装読んでみると常に有効桁数18あたりで最近接丸めした結果になるっぽい(真面目に読む気は無いので以下略)。

以下のようにdoubleの値の仮数部を前後に1bitづつずらした値をモード0で変換した結果を比較すると判りやすいだろうか。

$ cat >unko.c
#include <stdio.h>
#include <string.h>
#include "gdtoa.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;
 	char *head, *tail;
	int exp, neg;

	for (i = 0; i < arraycount(s); ++i) {
		memcpy(u.bvalue, s[i], sizeof(u.bvalue));
		head = dtoa(u.dvalue, 0, 0, &exp, &neg, &tail);
		printf("dvalue[%g], head:[%s]\n", u.dvalue, head);
	}
}
^D
$ make unko
$ ./unko
dvalue[0.1], head:[9999999999999999]
dvalue[0.1], head:[1]
dvalue[0.1], head:[10000000000000002]

まず%gの場合デフォルトの有効桁数6で切上げ発生するのでどれも0.1として表示されてるんだけど、モード0だとどれも違う値として変換されてる事がお判りいただけるかと *2

よって今回の%gの実装にこのモード0は使わず

  1. %eとおなじモード2を使う
  2. precが0の場合は1として扱う
  3. %eの場合の精度は小数点以下、%gの場合は整数部の1桁目も含むのでprecは-1ずれる
  4. 有効桁数までのゼロ埋めは不要(#まで実装するなら要るけど)

ちゅーかんじで実装するってわけよ。

コードは長くなるけど前々回の%f、前回の%eのコードを理解してればprecの調整とゼロ埋め以外は同じコードだと判るかと思うので、一気に全部貼る。

#include <assert.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#if defined(HAVE_MATH_H)
#include <math.h>
#endif

#include "gdtoa.h"

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

struct dtoatestcase {
	double dvalue;
	const char *svalue;
};

#define DTOATEST(value) { .dvalue = value, .svalue = #value }
#define debug	printf

#define PRECISION	6
#define MAXEXPSIZ	4

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
#if defined(HAVE_MATH_H)
		DTOATEST(INF),
		DTOATEST(NAN),
#else
		DTOATEST(1.0/0.0),
		DTOATEST(0.0/0.0),
#endif
		DTOATEST(0.0),
		DTOATEST(1.0),
		DTOATEST(12.0),
		DTOATEST(123.0),
		DTOATEST(10.0),
		DTOATEST(100.0),
		DTOATEST(1000.0),
                DTOATEST(1.1),
                DTOATEST(1.12),
                DTOATEST(1.123),
                DTOATEST(12.1),
                DTOATEST(12.12),
                DTOATEST(12.123),
                DTOATEST(123.1),
                DTOATEST(123.12),
                DTOATEST(123.123),
                DTOATEST(0.1),
                DTOATEST(0.01),
                DTOATEST(0.002),
		DTOATEST(1.123456789),
		DTOATEST(11.23456789),
		DTOATEST(112.3456789),
		DTOATEST(1123.456789),
		DTOATEST(11234.56789),
		DTOATEST(112345.6789),
		DTOATEST(1123456.789),
		DTOATEST(11234567.89),
		DTOATEST(112345678.9),
		DTOATEST(1123456789.0),
		DTOATEST(11234567890123.0),
	};
	size_t i;
	char *head, *tail, expstr[MAXEXPSIZ + 2];
	int prec, gprec, exp, neg, len, zeros, expsig;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec <= PRECISION; ++prec) {
			debug("dvalue:[%1$s], gprec:[%2$d], expected:[%3$.*2$g], result:[",
			    t[i].svalue, prec, t[i].dvalue);
			gprec = prec;
			if (gprec == 0)
				gprec = 1;
			head = dtoa(t[i].dvalue, 2, gprec, &exp, &neg, &tail);
			if (head == NULL)
				abort();
			len = tail - head;
			if (neg) {
				if (putc('-', fp) == EOF)
					abort();
			}
			if (exp == 9999) {
				if (fwrite(head, 1, len, fp) != len)
					abort();
			} else if (exp > -4 && exp <= gprec) {
				if (len <= exp) {
					if (fwrite(head, 1, len, fp) != len)
						abort();
					for (zeros = len; zeros < exp; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
				} else if (exp > 0) {
					if (fwrite(head, 1, exp, fp) != exp)
						abort();
	 				if (putc('.', fp) == EOF)
						abort();
					len -= exp;
                        		if (fwrite(&head[exp], 1, len, fp) != len)
                                		abort();
				} else {
 					if (fwrite("0.1", 1, 2, fp) != 2)
						abort();
					for (zeros = exp; zeros < 0; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
                        		if (fwrite(head, 1, len, fp) != len)
						abort();
				}
			} else {
				if (putc(*head, fp) == EOF)
					abort();
				if (--len > 0) {
					if (putc('.', fp) == EOF)
						abort();
					if (write(&head[1], 1, len, fp) != len)
						abort();
				}
				if (--exp >= 0) {
					assert(exp <= DBL_MAX_EXP);
					expsig = '+';
				} else {
					assert(exp >= DBL_MIN_EXP);
					expsig = '-';
					exp = -exp;
				}
				tail = &expstr[sizeof(expstr)];
				if (exp > 9) {
					do {
						*--tail = '0' + (exp % 10);
						exp /= 10;
					} while (exp > 9);
					*--tail = '0' + exp;
				} else {
					*--tail = '0' + exp;
					*--tail = '0';
				}
				*--tail = expsig;
				*--tail = 'e';
				len = &expstr[sizeof(expstr)] - tail;
				if (fwrite(tail, 1, len, fp) != len)
					abort();
			}
			freedtoa(head);
			debug("]\n");
		}
	}
}

もちろんコードを整理すれば%fや%gのコードと共有できるけどそれは各人で勝手にやってどうぞ。

@次回予告

露骨に手抜きになってきたけど、次回はいつになるか知らないけどhdtoaで%aを実装する予定。

*1:元天才塾にありそうな部活である。
*2:doubleの値を直接つくるには共用体を使うとよい、Nの場合にはmachine/ieee.hにieee_double_uそしてglibcでもmath_private.hにieee_double_shape_typeなんてのがご用意されてるはず。 なお今回はめんどくせえのでchar[8]との共用体で書いたのでENDIAN無視したコードなのでそこだけ注意な。