先週は平方根関数をライブラリとして実装する方法をひとつ紹介しました。その方法は、デフォルトの丸めモード（最近接丸め）では正しく動作しますが、方向丸めの場合にはうまくいかないことが判明しました。
その中で、逆数平方根をニュートン・ラフソン法で求める、invSqrt関数を実装しました。その実装は、丸め誤差に気を使っていない、ナイーブな実装になっています。それが原因で方向丸めの場合にうまくいかないのかを確かめるため、この関数の性質について調べます。
invSqrt関数の定義
以下で他の実装についても議論するので、invSqrt0とします。
float apploxInvSqrt( float x ) {
static_assert( sizeof(float) == sizeof(uint32_t) );
uint32_t i;
std::memcpy( &i, &x, sizeof i );
i = 0x5f3759df - ( i >> 1 );
float f;
std::memcpy( &f, &i, sizeof f );
return f;
}
float invSqrt0( float x ) {
// overflows if x <= 0x1.d1cc6p-129
assert( isnormal(x) && x > 0 );
float f = apploxInvSqrt(x);
f *= 1.5f - 0.5f * x * ( f * f );
f *= 1.5f - 0.5f * x * ( f * f );
f *= 1.5f - 0.5f * x * ( f * f );
return f;
}
逆数平方根をニュートン・ラフソン法で求める方法の性質
- （無限精度で計算した場合、）以下の数式で求める場合、収束は単調になる。
invSqrt0実装の性質
invSqrt0の実装方法では、係数1.5f - 0.5f * x * (f*f)の計算を無限精度でやった後一度だけ丸めたとしても、返り値に0.5ulpを超える誤差が乗ることがある。
たとえば、x = 0x1.00127cp+0fの場合、初期予測は0x1.eebaa8p-1f、一回加速後は0x1.ff1084p-1f、二回加速後は0x1.ffecf4p-1fとなりますが、この時つぎの係数は正確に計算すると0x1.00004882ac8fcp+0程度になります。これを丸めると0x1.000048p+0fになるが、これを乗じると0x1.ffed84p-1fになってしまいます。
切り上げで丸めると0x1.00004ap+0fを得ることができますが、これを乗じた場合は0x1.ffed88p-1fになってしまいます。つまり、真の値0x1.ffed850031aa1dp-1を丸めた0x1.ffed86p-1になるような係数は、そもそも単精度浮動小数点数の中に存在しません。
これは、求めるべき値の仮数部が2に近く、精度が必要であるのに、係数の仮数部は1に近く、精度が一桁分足りなくなってしまっていることが原因です。
invSqrt0の実装方法では、複数回の丸めが原因で、返り値に1.9ulpを超えるような誤差が乗ることがある。
たとえば、x = 0x1.08fd12p+0fの場合、二回加速後は0x1.f73d9ap-1fになります。この時、真の係数は0x1.00001a813f05dp+0程度です。計算過程を追ってみると、f2*f2 == 0x1.eea19p-1fですが、これは0x1.eea190ff6052p-1を丸めたものであり、0.5ulp近く切り下げられてしまっています。さらに0.5f * x * (f2*f2) == 0x1.ffff94p-2fですが、これは0x1.ffff95fb03e8bp-2を丸めたものであり、再び0.5ulp近く切り下げられてしまいました。
次に、1.5f - 0x1.ffff94p-2fが計算され、結果は0x1.00001bp+0となるはずですが、偶数丸めが発動して0.5ulp切りあがります。複数回丸めてはいけないという典型例になっています。
全ての丸めが同じ方向に働いたため、誤差が大きくなってしまいました。 ちなみに、この例では、誤差は1.9485ulpに達しています。
- 丸め誤差が原因で、何度加速しても収束しないことがある。
f *= 1.5f - 0.5f * x * ( f * f );
の加速を複数回繰り返しても、誤差が原因で二つの値を振動することがあります。
丸め誤差を減らすために最後にstd::fmaを使うinvSqrt1実装の性質
float kasoku( float x, float y ) {
float hh = 0.5f - (0.5f*x) * (y * y);
float ret = std::fma(y, hh, y);
return ret;
}
float InvSqrt1( float x ) {
// overflows if x <= 0x1.d1cc6p-129
assert( isnormal(x) && x > 0 );
float f = apploxInvSqrt(x);
f *= 1.5f - 0.5f * x * ( f * f );
f *= 1.5f - 0.5f * x * ( f * f );
f = kasoku( x, f );
return f;
}
このコードは、逆数平方根の収束公式と丸め誤差のconv_2nd_B関数を参考にしました。
最終回の加速の計算で
std::fmaを一回使う方法を用いても、誤差は1ulpには収まらない。上記のコードを用いてもう一度加速すると却って誤差が増える。
上記のコードを用いてさらにもう一度加速すると、誤差は1ulpに収まる*1。
丸め誤差を取り除くために追加で二回も加速しないといけないようでは、もはや加速とは呼べないような気がします。
二つのfloatを用いて極力丸め誤差を減らすinvSqrt2実装の性質
float kasoku( float x, float y ) {
float y2 = y*y;
float yc = std::fma( y, y, -y2 );
float h = std::fma( -0.5f*x, y2, 0.5f );
float hc = std::fma( -0.5f*x, yc, h );
float ret = std::fma( y, hc, y );
return ret;
}
float InvSqrt2( float x ) {
// overflows if x <= 0x1.d1cc6p-129
assert( isnormal(x) && x > 0 );
float f = apploxInvSqrt(x);
f *= 1.5f - 0.5f * x * ( f * f );
f *= 1.5f - 0.5f * x * ( f * f );
f = kasoku( x, f );
return f;
}
このコードは、逆数平方根の収束公式と丸め誤差のconv_2nd_C関数を参考にした。二つのfloatを用いて48bit分の精度を確保する試みです。
この極力丸め誤差を減らした加速方法を三回目の加速に用いても、必ずしも真値を丸めた値になるとは限らない。ただし、誤差は0.5004ulpを超えない*2。
それどころか、上記のコードを用いた場合でも、何度加速しても収束しない場合がある。
x = 0x1.13e070p+1fの場合、真の逆数平方根は0x1.5cc0a9000000bp-1程度であり、48bit精度ではどちらに丸めてよいか判別がつきません。実際、0x1.5cc0a8p-1fと0x1.5cc0aap-1の間で振動します。
まとめ
調べてみましたが、あまり得るものがありませんでした……。最後に出てきた、0x1.13e070p+1fの真の逆数平方根が0x1.5cc0a9000000bp-1であるという例はかなり凶悪ですね。
逆数平方根を正確に求めようとすると結構泥沼のようです。
mySqrt関数の改善
もっとも単純なinvSqrt0の方法でも誤差は2ulp以内であり、std::fesetround(FE_DOWNWARD);されている時に先週公開したmySqrt関数をこれを使っても、ごく一部の平方数に対する平方根をぴったり1ulp小さく算出してしまう以外の間違いは発生しないようです。
invSqrt1を用いるとそのような間違いの頻度はさらに下がります。invSqrt2を用いた場合、invSqrt0よりは減りますがinvSqrt1よりは増えるという謎現象が発生します。
