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

2021/04/09(Fri)

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

(追記) 第二引数についてとんでもない勘違いしてたのと、それ以降の引数の番号がズレてたので修正しました。

数学は高二で脱落した「さんすうが苦手か…!?」な人だしCPUとFPUとGPUそしてxPUの違いもぜんぜんわからない俺は雰囲気でコードを書いてるマンなので、浮動小数数数数点点点点数数の扱いでIEEEEEEEE754の仕様を読もうとするだけで猛烈な頭痛に襲われ寝込むので、そもそもプログラマ名乗ってたこと自体間違いなのである。

なおワイのように不動小数と書いたりIEEEのEの数がいくつだったも数えられないのでいつも大目に書いてる莫迦であっても、倍賞美津子度浮動小数点数計算ハコテン操作できるライブラリが表題のgdtoaである、それロン!(脱衣シーン鑑賞中)

ライセンスがおおよそ2-clause BSDLと同等で緩いから利用に問題は無いし、それに自分でIEEE754の実装して楽しいと思える人類は世界で3人くらいなので、世のほとんどのlibc実装やLL言語もこのgdtoa使って実装されてるのだが、ソースがOpenSSLすら読みやすく思えてくる古代のインタフェースも洗練されてないので、いつか捨ててやりてえと思ってるが多分死ぬまで無理そう。

このgdtoaってのは元々dtoa.cというdoubleを文字列に変換するCソース単体で配布されてたので、まずはそいつの使い方。

#include <stddef.h>
#include <stdio.h>
#include "gdtoa.h"
int
main(void)
{
	double dvalue = 1.123456789;
	char *head, *tail;
	int exp, neg;

	head = dtoa(dvalue, 2, 8, &exp, &neg, &tail);
	if (head == NULL)
		abort();
	if (exp == 9999)
		...
	...
	freedtoa(head);
}

このdtoa()という関数を呼ぶだけで文字列に変換してくれる。

まず二番目で指定するのは変換モードですな。

        mode:
                0 ==> shortest string that yields d when read in
                        and rounded to nearest.
                1 ==> like 0, but with Steele & White stopping rule;
                        e.g. with IEEE P754 arithmetic , mode 0 gives
                        1e23 whereas mode 1 gives 9.999999999999999e22.
                2 ==> max(1,ndigits) significant digits.  This gives a
                        return value similar to that of ecvt, except
                        that trailing zeros are suppressed.
                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.
                4,5 ==> similar to 2 and 3, respectively, but (in
                        round-nearest mode) with the tests of mode 0 to
                        possibly return a shorter string that rounds to d.
                        With IEEE arithmetic and compilation with
                        -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
                        as modes 2 and 3 when FLT_ROUNDS != 1.
                6-9 ==> Debugging modes similar to mode - 4:  don't try
                        fast floating-point estimate (if applicable).

                Values of mode other than 0-9 are treated as mode 0.

                Sufficient space is allocated to the return value
                to hold the suppressed trailing zeros.

ざっと簡単に説明すると

ということで普通に2と3の使い方だけ知ってればほぼ事足りると思われる(おい。

誤差の丸め方についてはヘッダファイルに定義があるけど

enum {  /* FPI.rounding values: same as FLT_ROUNDS */
        FPI_Round_zero = 0,
        FPI_Round_near = 1,
        FPI_Round_up = 2,
        FPI_Round_down = 3
        };

の4つ、これはIEEE754-1985で定義されてるもので

いま日本語訳として適切なもの探そうとして、普通に「マルコメ」って検索欄に打ち込んじゃって味噌汁出てきたし、出てきた用語も最近セッとか性の無限大と空目するし、ついには「丸め」が「攻め」に見えだしてもうダメだ。 なんせ∞という記号をみるだけで脳内でヨコハチ無限大と読んでしまうし三浅一深とか秘技卍ハーケンクロスとか以下略。

ちなみにIEEE754-2008で追加された

は未対応なはず、なんせ最終更新が1998年だからな…どっかに対応版とかあったら教えて。

そんで三番目の引数は有効桁数の指定、printf(3)なんかだと

$ cat >unko.c
#include <stdio.h>
int
main(void)
{
	printf("%.8f\n", 1.123456789);
}
^D

などと精度として指定するものと思えばよろし、ただしこれを実行すると

$ make unko
$ ./unko
1.12345679

と、有効桁数以降を丸めた結果が表示されるのだけど、dtoaが返す文字列はこのままでは無いことに注意(後述)。

四番目の引数は 山椒渡しポインタ渡しになってるけど、これはdtoaからの 返り血戻り値で、exponentつまり指数が入ってくる。 ワイは3本早く人間になりたい~詳しいことはこれも後述。 ちなみに9999が返された場合はdoubleの値がINF(Infinity)あるいはNAN(Not A Number)だったという意味を持ち、変換結果は

のどっちかになる。なおmath.hにあるisinf(3)は主にVAXとかVAXとかVAX方面がサポートしておらず移植性を欠くので未だに現役の某OSで困るのだが、gdtoaではどうなってるかというと

#if defined(IEEE_Arith) + defined(VAX)
#ifdef IEEE_Arith
        if ((word0(&d) & Exp_mask) == Exp_mask)
#else
        if (word0(&d)  == 0x8000)
#endif
                {
                /* Infinity or NaN */
                *decpt = 9999;
#ifdef IEEE_Arith
                if (!word1(&d) && !(word0(&d) & 0xfffff))
                        return nrv_alloc("Infinity", rve, 8);
#endif
                return nrv_alloc("NaN", rve, 3);
                }
#endif

この「#if defined(IEEE_Arith) + defined(VAX)」という謎のifdefで非常に混乱したのだがVAXとIEEE_Arithが同時にdefinedになるはずが無いのでやっぱりInfinityはサポートしてないちゅーこと。

ここではたと思い出したのだけど、昔N/i386のstrtold(gdtoaのstrtopx.c)が正しくINFを変換してくれないという記事を 書いて、まぁ俺の専門分野じゃねーからとsend-prだけしといたんだが、もしかしてgdtoaが正しくてN/i386のsrc/lib/libc/arch/i386/gen/infinityl.cの方が間違ってるような気がしてきたゾ。 とはいえ放置プレイ喰らってそのまま埋もれてるし、俺もいまさらi386なんて使ってないからどうでもいいのだけどね…

よく考えたらx86_64も同じだったわ、ローカルパッチ外したら

$ cat >unko.c
#include <math.h>
#include <stdlib.h>
#include <assert.h>

int
main(void)
{
        assert(isinf(strtold("INF", NULL)));
}
^D
$ make unko
$ ./unko
assertion "isinf(strtold("INF", NULL))" failed: file "unko.c", line 8, function "main"
Abort (core dumped)

とまぁ盛大に死んでくれよる。

オレオレN6ではさっきの記事の差分をstrtopx.cに適用してるのだけど、もしかしてNのINF定義(src/lib/libc/arch/*/gen/infinityl.c)

[i386]
const union __long_double_u __infinityl =
        { { 0, 0, 0, 0, 0, 0, 0, 0x80, 0xff, 0x7f, 0, 0 } };

[x86_64]
const union __long_double_u __infinityl =
        { { 0, 0, 0, 0, 0, 0, 0, 0x80, 0xff, 0x7f, 0, 0, 0, 0, 0, 0 } };

の方が間違ってるとしたらアレでアレ…あとでまじめに仕様調べて正しい修正を考えようか…

話を戻して五番目の引数、こいつも察しのとおり戻り値で、doubleの値が

が返される、よく負のオーラ滲み出てるといわれます天井のシミの顔がいつも攻撃してくるんです。

六番目の引数以下同文、dtoaの戻り値である文字列の末尾を指すポインタが返ってくる。 この仕様だけ見ると

的な事を想像してしまうのだけど、コード読む限りではちゃんとヌルヌルターミネーターされてるんだよな。

        if (s == s0) {                  /* don't return empty string */
                *s++ = '0';
                k = 0;
        }
        *s = 0;
        *decpt = k + 1;
        if (rve)
                *rve = s;
        return s0;
        }

なので使わなくてもいい、まぁケツの方から処理したい時に使えばいいんじゃないでしょうか、剃刀とか。

このdtoaの戻り値の文字列は必ずfree(3)ではなくfreedtoaで開放する必要がある、というのもgdtoa元々はdtoa_resultというグローバル変数を使い回すのでスレッドアンセーフだったのだ。

#ifndef MULTIPLE_THREADS
 char *dtoa_result;
#endif

 char *
#ifdef KR_headers
rv_alloc(i) size_t i;
#else
rv_alloc(size_t i)
#endif
{
        size_t j;
        int k, *r;

        j = sizeof(ULong);
        for(k = 0;
                sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
                j <<= 1)
                        k++;
        r = (int*)(void*)Balloc(k);
        if (r == NULL)
                return NULL;
        *r = k;
        return
#ifndef MULTIPLE_THREADS
        dtoa_result =
#endif
                (char *)(void *)(r+1);
        }

なもんで

という対応を入れよったからなのだ、うーんこの。ところでfree(3)のラッパーって流しのMCかなんかか?

よって安全側に倒すために

	size_t len;
	char *head, *tail, *tmp;
	int exp, neg;

	head = dtoa(dvalue, 2, 8, &exp, &neg, &tail);
	if (head == NULL)
		abort();
	len = tail - head;
	tmp = strdup(head);
	if (tmp == NULL)
		abort();
	freedtoa(head);
	head = tmp;
	tail = tmp + len;

とでもコピーしちまった方がいいかもね、こういうところがパラノイアっていわれるねん。

@次回予告

以前printf(3)のウン国際化対応でさんざんgdtoa呼び出してるスパゲッティーを超えた鬱コードを触ってもう二度と触りたくないと思ってたのだが、暇潰しがてらのJSON parser for Cの実装でprintf("%g")をサポートしない環境を考えたら結局またgdtoa触る羽目になってしまった。 次は戻り値の文字列がどのような形式で返されるのか、そして指数と正負をつかってどのように普段目にする%fとか%eそして%gまたまた%aの形式に変換するのかを説明する予定、いつになるかは知らん。

[プログラミング] gdtoaのstrtold(strtopx.c)のバグはすでに本家で修正済であった

さっきの記事、最新版のgdtoaみたらワイの当初の修正で正しかったようでgdtoa側に同じ修正が入ってる事が判明、というか更新があったとかしらそん。そんでNもマージ済なのでN7以降では直ってるというちゅことか。 それでもワイのsend-prから5年くらいは放置されたわけでまぁ世の中isinf(3)が動かなくてもいかに何の問題も無いかということだ、というかgdtoaの作者さんgithubで管理してくだち!県北の土手でプルリク盛りあおうぜ。