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

2021/04/11(Sun)

[プログラミング] gdtoaの使い方(その2)

前回の続き、あくまで使い方であってgdtoa(のソースコードの)歩き方ではない…はず…

本当は先に用語の説明を書こうと思ってたのだが、文章力がボブ・ディラン並みなので後回しにしてコード例を先に置いておく。

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

まずはgdtoaのパラメーターである

  • 第二引数の変換モード
  • 第三引数の有効桁数

に何をセットするかなんだけど、モードについては *1前回も引用したコメント

                3 ==> through ndigits past the decimal point.  This
                        gives a return value similar to that from fcvt,
                        except that trailing zeros are suppressed, and
                        ndigits can be negative.

いつもの超訳だと

	3 ==>	小数点以降を有効桁数とする、これは%f変換とよく似た結果となる。
		ただし末尾の連続するゼロは取り除かれるし、有効桁数は負の値になることもある。

と書かれているので3を使うことにする。有効桁数は負の値ってのが意味不明かもだが後で説明でして?そいつはコトだ。

そんでCの仕様だと精度で明示的に指定しなかった場合は小数点以下の有効桁数は6となる

If the precision is missing, it shall be taken as 6;

よって第三引数はそれで。

#include "gdtoa.h"

#define PRECISION	6

int
main(int argc, char *argv[])
{
	double dvalue;
	char *head, *tail;
	int prec = PRECISION, exp, neg, len;

	head = dtoa(dvalue, 3, prec, &exp, &neg, &tail);
	len = tail - head;
	...
}

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

ここ読んでるすうがくやってる理系のコンピーターおじいちゃんたち(なお想定読者層はシティポップ好きゆるふわカメラ女子)は有効数字(Significant Digit)知ってると思うけど、dtoaはそっからさらに右側のゼロを削った形式を返す。 なんか他にこれ呼び名あるんだっけ、数学は理系クラス男ばっかりだから高2で投げ出したからな。 まぁ要は左や右の旦那さま(政治)とゼロは全部消しちまったぜってことです。

そんで削られた左側のゼロおよび小数点はexpから復元、そして右側のゼロはlenおよびprecから復元できるという寸法。

たとえば「0.00123456789000の有効数字」というと、左側のゼロ及び小数点を無視し「1234567890000」の12桁となるけど、gdtoaではさらに右側も無視して「123456789」の9桁になる。 さらに今回のサンプルコードでは「有効桁数は小数点以下6桁」としたので、7桁目で切上げを行うので最終的には「1235」が返ってくる文字列となる。

そんで小数点の位置から数えて文字列の先頭は

0
.
0 <= [ 0]
0 <= [-1]
1 <= [-2]
2
3
5

と表現できるけどこれがexpとして返される値となる。マイナスになるってのがさっきのコメントにあった「有効桁数が負の値」の意味だと思われる、作者はえいごがにがてか…!?

じゃあ実際に右や左のゼロと小数点を復元して普段目にする形式に変換するロジックを書いていこうか。

@A. expの値が9999となるケース

これは前回も軽く触れたけど、値がINFあるいはNANだった場合のケース。

dvalue len exp neg head
INF 8 9999 0 "Infinity"
NAN 3 9999 0 "NaN"
-INF 8 9999 1 "Infinity"
-NAN 3 9999 1 "NaN"

このケースはそのままheadを出力すればいい。

#include <assert.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

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
#if defined(HAVE_MATH_H)
		DTOATEST(INF),
		DTOATEST(NAN),
		DTOATEST(-INF),
		DTOATEST(-NAN),
#else
		DTOATEST(1.0/0.0),
		DTOATEST(0.0/0.0),
		DTOATEST(-1.0/0.0),
		DTOATEST(-0.0/0.0),
#endif
	};
	size_t i;
	char *head, *tail;
	int prec, exp, neg, len;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec < PRECISION; ++prec) {
			debug("dvalue:[%s], prec:[%d], expected:[%.*f], result:[",
			    t[i].svalue, prec, prec, t[i].dvalue);
			head = dtoa(t[i].dvalue, 3, prec, &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) != 3)
					abort();
			}
			freedtoa(head);
			debug("]\n");
		}
	}
}

もしmath.hが存在しない環境であれば 過去記事にも書いたように

  • INF … (1.0/0.0)
  • NAN … (0.0/0.0)

で作れる(お前がVAXでも使ってない限りは)、そいやVisual C/C++だとこれ定数で書けない問題って直ったのだろうか。

なお仕様ではINFおよびNANは

The F conversion specifier produces "INF", "INFINITY", or "NAN" instead of "inf", "infinity", or "nan", respectively.

と変換しロッテ書かれてるので、upper/lower混在となるgdtoaの出力は正しくないことに注意。

あと今Nのprintf(3)が負のNANを正として出力するバグを発見したような気がするが気にしない気にしない。

@B. len == expとなるケース

判りやすいように例を挙げると以下のケース。

dvalue len exp neg head
0.0 1 1 0 "0"
1.0 1 1 0 "1"
12.0 2 2 0 "12"
123.0 3 3 0 "123"
-0.0 1 1 1 "0"
-1.0 1 1 1 "1"
-12.0 2 2 1 "12"
-123.0 3 3 1 "123"

最後のやつを例にとると、小数点の位置から数えて文字列の先頭は

-
1 <= [3]
2 <= [2]
3 <= [1]
.
0

と表現されるのでexp = 3が返ってくるわけすわ。

このケースでは

  • 小数点以下は全てゼロ(フレドリック・ブラウン著)
  • よってINF/NAN同様にそのままheadは全部出力してしまってよろしい
  • あとは小数点以下を有効桁数までゼロで埋めればいい

というコードになる。

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

#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

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
		DTOATEST(0.0),
		DTOATEST(1.0),
		DTOATEST(12.0),
		DTOATEST(123.0),
		DTOATEST(-0.0),
		DTOATEST(-1.0),
		DTOATEST(-12.0),
		DTOATEST(-123.0),
	};
	size_t i;
	char *head, *tail;
	int prec, exp, neg, len, zeros;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec <= PRECISION; ++prec) {
			debug("dvalue:[%s], prec:[%d], expected:[%.*f], result:[",
			    t[i].svalue, prec, prec, t[i].dvalue);
			head = dtoa(t[i].dvalue, 3, prec, &exp, &neg, &tail);
			if (head == NULL)
				abort();
			len = tail - head;
			if (neg) {
				if (putc('-', fp) == EOF)
					abort();
			}
			if (len == exp) {
				if (fwrite(head, 1, len, fp) != len)
					abort();
				if (prec > 0) {
					if (putc('.', fp) == EOF)
						abort();
					for (zeros = 0; zeros < prec; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
				}
			}
			freedtoa(head);
			debug("]\n");
		}
	}
}

上のコード例のとおり、IEEE 754には負のゼロ(-0.0)があるのだけど、よくある間違いとして

	double dvalue = -0;

と書いたりすると負のゼロが存在しえないintのリテラルからのキャストとなるのでお前は書いたはずの符号を失って死ぬ。

そして結果をみてもうお気づきでしょうが、123.0と-123.0でdtoaの戻り値が違うのはnegが0か1かだけであとは全く一緒。

というのもdoubleの内部表現は

  • 符号部(sign) … 1bit
  • 仮数部(fraction) … 52bit
  • 指数部(exponent) … 11bit

であるので、この1bitの違いでしかないので当然っちゃー当然、なので以降のケースでは負の値は省略させてもらう。

そりゃ正か負かなんて符号があるかないだけだもんな、もし銀行預金額がdoubleで記録されてたら記憶媒体の1bitが反転しただけで借金が貯金に反転するなんて夢あるよね…やはりECCメモリの採用は悪…化けろもっと化けろ…

まぁ金融系はふつー10進数使ってて2進数なんぞ使わんはずなので以下略、これは別に1bit化けるだけ云々という話ではなく誤差の問題でだぞ。

@C. len < expとなるケース

こいつも例を挙げると

dvalue len exp neg head
10.0 1 2 0 "1"
100.0 1 3 0 "1"
1000.0 1 4 0 "1"

最後のやつを例にとると、小数点の位置から数えて文字列の先頭は

1 <= [4]
0 <= [3]
0 <= [2]
0 <= [1]
.
0

と表現されるのでexp = 4が返る。

このケースでは

  • 小数点の後だけでなく前にも補完が必要なゼロが存在する
  • よってheadをすべて出力した後0の位まで(つまりexp - len個の)ゼロ埋めする必要がある
  • 小数点以下を有効桁数precまでゼロで埋めるのは一緒

ってのがlen == expのケースとの違いですな。

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

#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

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
		DTOATEST(10.0),
		DTOATEST(100.0),
		DTOATEST(1000.0),
	};
	size_t i, len;
	char *head, *tail;
	int prec, exp, neg, zeros;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec <= PRECISION; ++prec) {
			debug("dvalue:[%s], prec:[%d], expected:[%.*f], result:[",
			    t[i].svalue, prec, prec, t[i].dvalue);
			head = dtoa(t[i].dvalue, 3, prec, &exp, &neg, &tail);
			if (head == NULL)
				abort();
			len = tail - head;
			if (neg) {
				if (putc('-', fp) == EOF)
					abort();
			}
			if (len < exp) {
				if (fwrite(head, 1, len, fp) != len)
					abort();
				for (zeros = len; zeros < exp; ++zeros) {
					if (putc('0', fp) == EOF)
						abort();
				}
				if (prec > 0) {
					if (putc('.', fp) == EOF)
						abort();
					for (zeros = 0; zeros < prec; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
				}
			}
			freedtoa(head);
			debug("]\n");
		}
	}
}

ここで勘のいい人は気づいただろうけど

				for (zeros = len; zeros < exp; ++zeros) {
					if (putc('0', fp) == EOF)
						abort();
				}

このループはlen == expの場合には回らない、よってlen == expとlen < expのコードはほぼ同じなのだ。

だからコードを短くしたければ

			if (len == exp) {
				block B
			} else if (len < exp) {
				block C
			}

の条件分岐は

			if (len <= exp) {
				block C
			}

としておまとめてしまってよろしいと思います。

@D. len > exp かつ exp > 0のケース

はいはい、サンプルね。

dvalue len exp neg head
1.1 2 1 0 "11"
1.12 3 1 0 "112"
1.123 4 1 0 "1123"
12.1 3 2 0 "121"
12.12 4 2 0 "1212"
12.123 4 2 0 "12123"
123.1 4 3 0 "1231"
123.12 5 3 0 "12312"
123.123 6 3 0 "123123"

最後のやつを例にとると、小数点の位置から数えて文字列の先頭は

1 <= [3]
2 <= [2]
3 <= [1]
.
1
2
3

と表現されるのでexp = 3が返る、もう覚えたね?簡単だね!

このケースでは

  • headを小数点の前つまり長さexpだけ出力する
  • headの残りlen - exp分は小数点出力後に
  • precまでゼロ埋め

というロジックでおk。

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

#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

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
                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),
	};
	size_t i;
	char *head, *tail;
	int prec, exp, neg, len, zeros;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec <= PRECISION; ++prec) {
			debug("dvalue:[%s], prec:[%d], expected:[%.*f], result:[",
			    t[i].svalue, prec, prec, t[i].dvalue);
			head = dtoa(t[i].dvalue, 3, prec, &exp, &neg, &tail);
			if (head == NULL)
				abort();
			len = tail - head;
			if (neg) {
				if (putc('-', fp) == EOF)
					abort();
			}
			if (len > exp && exp > 0) {
				if (fwrite(head, 1, exp, fp) != exp)
					abort();
				if (prec > 0) {
	 				if (putc('.', fp) == EOF)
						abort();
					len -= exp;
                        		if (fwrite(&head[exp], 1, len, fp) != len)
                                		abort();
					for (zeros = len; zeros < prec; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
				}
			}
			freedtoa(head);
			debug("]\n");
		}
	}
}

@E. len > expかつexpが0以下のケース

dvalue len exp neg head
0.1 1 0 0 "1"
0.01 1 -1 0 "1"
0.001 1 -2 0 "1"

最後のやつを例にとると、小数点の位置から数えて文字列の先頭は

0
.
0 <= [ 0]
0 <= [-1]
1 <= [-2]

と表現されるのでexp = -2が返る、もう完璧!

よってロジックは

  • 先頭のゼロと小数点を出力
  • -exp個ゼロ埋めをする
  • headを全て出力する
  • precまで残りゼロ埋め
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#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

int
main(int argc, char *argv[])
{
	static const struct dtoatestcase t[] = {
                DTOATEST(0.1),
                DTOATEST(0.01),
                DTOATEST(0.002),
	};
	size_t i;
	char *head, *tail;
	int prec, exp, neg, len, zeros;
	FILE *fp = stdout;

	for (i = 0; i < arraycount(t); ++i) {
		for (prec = 0; prec <= PRECISION; ++prec) {
			debug("dvalue:[%s], prec:[%d], expected:[%.*f], result:[",
			    t[i].svalue, prec, prec, t[i].dvalue);
			head = dtoa(t[i].dvalue, 3, prec, &exp, &neg, &tail);
			if (head == NULL)
				abort();
			len = tail - head;
			if (neg) {
				if (putc('-', fp) == EOF)
					abort();
			}
			if (len > exp && exp <= 0) {
 				if (putc('0', fp) == EOF)
					abort();
				if (prec > 0) {
 					if (putc('.', fp) == EOF)
						abort();
					for (zeros = exp; zeros < 0; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
                        		if (fwrite(head, 1, len, fp) != len)
						abort();
					for (zeros = len - exp; zeros < prec; ++zeros) {
						if (putc('0', fp) == EOF)
							abort();
					}
				}
			}
		}
		freedtoa(head);
		debug("]\n");
	}
}

@次回

はいこれでprintf(3)の%fとほぼ同等の変換がgdtoaで実現できるようになりました、まぁ他にも書式指定子による幅指定とか左寄せゼロorスペース埋めとかも考慮する必要があるのだけど、そこいらへんは割愛します。

次回は%eの変換について説明する予定、いつになるかは知らん。

ただし最もめんどくせえPOSIX localeによる国際化機能すなわちLC_NUMERICにおける

  • decimal_point … 小数点記号
  • thousands_sep … 桁区切り記号
  • grouping … 桁区切り記号の挿入位置のルール

についてもおそらくprintf(3)を自分で実装する機会でもなけりゃまず使い方とか知ることが無いだろうから、これについては機会があれば説明書きたいとは思っている。

ちなみにワイまだ幼児の頃にみた古いIBMのCMがお気に入りだったのだが、そのせいで日本人がアラビア数字と1000毎に桁区切り入れることに非常に違和感があるのだ。

なのでその結果さんすうが苦手にと思われる、おのれIBM!

*1:前回の説明でなぜか第二引数丸めモードの指定と勘違いして記事書いてたわ(訂正済)、わからないやっぱり俺は雰囲気でgdtoaを操作している。