hack のためのネタ帳, etc,,,

状況

症状としては以下のページにあるように、fenv.h の fesetround() で浮動小数点演算の丸めモードを変更しても、math.h 関連の関数に対して効果が得られない。

上記で、例示されている math.h の sqrt() は Cygwin の場合、以下のパッケージで提供される /bin/cygwin1.dll 内に実装されている。 こいつの source は以下のは、newlib-cygwin が source となっている。
現在のバージョンは cygwin-3.1.7-1

ここで、sqrt() は
  • cygwin-3.1.7-1.src/newlib-cygwin/newlib/libm/math/w_sqrt.c
で定義されてて、これは以下のファイルで定義されている double __ieee754_sqrt(double x) への wrapper である。
  • cygwin-3.1.7-1.src/newlib-cygwin/newlib/libm/math/e_sqrt.c

これを、拾ってきて
fesetround(roundmode);
x = sqrt(a);
y = __ieee754_sqrt(a);
のようなコードに対して
gcc myprog.c e_sqrt.c
みたいにしてコンパイルしてやると最適化の設定によって症状にブレが生じる。
cygwin1.dll に組み込まれている sqrt() は最適化が無効の場合丸めモードが機能せず、最適化を有効にすると丸めモードが機能する。
対して、自分でコンパイルした、__ieee754_sqrt() は、最適化が無効の場合、丸めモードが機能し、最適化を有効にすると丸めモードが機能しない。

cygwin1.dll は -O2 で最適化されているので、自分のプログラムを最適化無効で作成すると、これが呼ばれるため、丸めモードが機能しない。
一方、最適化有効にすると、stdc の math.h の関数は最適化により、cygwin1.dll のルーチンではなくもっと速いルーチンに差し替えられるらしく、こっちは fesetround() の効果が得られているという状況らしい。

原因

newlib が、丸め処理を自前で判定しているのが原因。
例えば、sqrt() の場合 89 行目
static	const double	one	= 1.0, tiny=1.0e-300;
とかいう定数を用いて、174-186行目
    /* use floating add to find out rounding direction */
	if((ix0|ix1)!=0) {
	    z = one-tiny; /* trigger inexact flag */
	    if (z>=one) {
	        z = one+tiny;
	        if (q1==(__uint32_t)0xffffffff) { q1=0; q += 1;}
		else if (z>one) {
		    if (q1==(__uint32_t)0xfffffffe) q+=1;
		    q1+=2; 
		} else
	            q1 += (q1&1);
	    }
	}
みたいにして丸めモードを判別して数値を補正している。
ここで、one-tiny とか one+tiny は加減算の左右の項が定数であるため、最適化が行われるとコンパイル時に値を計算して条件分岐を取り除いてしまうため、
実行時に fesetround() をいじっても効果が得られない。

なお、

みたいに言われているが、見ての通り Copyright (C) 1993 by Sun Microsystems, Inc 由来の古式ゆかしくも由緒正しいコードであり、計算方法の詳細についてもコメントで詳しく解説されてる。
アルゴリズム自体には問題はなくて、最適化の副作用でコケている状況だろう。

修正方法

tiny に valatile 型修辞子を付けて、
static	const volatile double	tiny=1.0e-300;
してやると、コンパイル時に計算されず、実行時に計算が行われるようになり、この問題を回避できた。

なお sqrt() 以外にも以下の箇所に tiny があったので、これらのファイルは影響を受けるはず(少なくとも atan2 は影響受けてた)。
$ tar xf newlib-cygwin-3.1.7.tar.bz2
$ grep -ER "tiny *=" . -n
./newlib-cygwin/newlib/libm/common/sf_scalbln.c:26:tiny   = 1.0e-30;
./newlib-cygwin/newlib/libm/common/sf_scalbn.c:34:tiny   = 1.0e-30;
./newlib-cygwin/newlib/libm/common/s_scalbln.c:32:tiny   = 1.0e-300;
./newlib-cygwin/newlib/libm/common/s_scalbn.c:76:tiny   = 1.0e-300;
./newlib-cygwin/newlib/libm/math/ef_atan2.c:24:tiny  = 1.0e-30,
./newlib-cygwin/newlib/libm/math/ef_sqrt.c:19:static    const float     one     = 1.0, tiny=1.0e-30;
./newlib-cygwin/newlib/libm/math/ef_sqrt.c:21:static    float   one     = 1.0, tiny=1.0e-30;
./newlib-cygwin/newlib/libm/math/e_atan2.c:51:tiny  = 1.0e-300,
./newlib-cygwin/newlib/libm/math/e_sqrt.c:89:static     const double    one     = 1.0, tiny=1.0e-300;
./newlib-cygwin/newlib/libm/math/e_sqrt.c:91:static     double  one     = 1.0, tiny=1.0e-300;
./newlib-cygwin/newlib/libm/math/sf_tanh.c:19:static const float one=1.0, two=2.0, tiny = 1.0e-30;
./newlib-cygwin/newlib/libm/math/sf_tanh.c:21:static float one=1.0, two=2.0, tiny = 1.0e-30;
./newlib-cygwin/newlib/libm/math/s_tanh.c:76:static const double one=1.0, two=2.0, tiny = 1.0e-300;
./newlib-cygwin/newlib/libm/math/s_tanh.c:78:static double one=1.0, two=2.0, tiny = 1.0e-300;
tiny を±する相手は one だったり pi だったりするので、相手側も volatile しておくべきかどうかは悩む所。

これ、少なくとも 5 年以上も放置されてるわけなんだが、ひょっとして Unit Test 組まれてないのだろうか?
報告は mailing list へと書いてあるので敷居が高いな。
英語はともかくとして、本垢で報告すると迷惑メール増えるしさ。

コメントをかく


「http://」を含む投稿は禁止されています。

利用規約をご確認のうえご記入下さい

Wiki内検索

フリーエリア

編集にはIDが必要です