IEEE浮動小数点数における演算では、丸め誤差が不可避です。特に、複数回の演算を繰り返すと丸め誤差が積もっていき、正確な値と大きく離れた答えを得てしまうことがあります。しかし、次の演算については、（数学的に）正確な値を求めた後、一回だけの丸めが発生することが、IEEE標準で規定されています。
浮動小数点数演算のできるCPUであれば、普通、四則演算や積和演算を行う命令を持っています。 しかし、平方根を正確に計算する命令を持たない命令セットも存在します。
そのような場合、平方根関数はライブラリ実装となるわけですが、どのように実装すれば要求を満たせるのでしょうか？
C++のstd::sqrtは正確に計算しているのか？
結論
しています。
標準の丸めモード、つまり最近接丸め（ぴったり半分ならば偶数へ丸める）を用いた時、単精度浮動小数点数floatについては全ての値で正しく丸められていることを確認しました。
確認は、Wandboxのclang7.0.0で行いました。
その他の性質
std::fmaとの関係
また、任意の単精度浮動小数点数xに対して、次のことが成り立つようです。
std::sqrt(x)は、std::abs(std::fma(s,s,-x))を最も小さくする正のsである。
これは自明ではありません。丸め誤差が0.499999999ulpであるほうに丸めるより、-0.500000001ulpであるほうに丸める方がstd::abs(std::fma(s,s,-x))が小さくなる、という場合があり得るわけです。そのようなぎりぎりの位置に平方根が存在することがないことは証明できるのでしょうか。四則演算の場合はそんな感じの定理があったような……？
これは実験的に確かめたことであり、少し怪しいです。特に、xが非正規化数である場合はそのようなsが複数個あることがあります（s*s-xを無限精度で計算できる場合には、sが正かつfloatで表現可能、という条件のもとそのようなsは一意に定まるはずですが、std::fma(s,s,x)がほぼ0（絶対値が最小の正の非正規化数よりも小さく、丸めると0.0fか-0.0fになる値）になり、精度が不足することが原因のようです）。
丸め誤差の最大値
丸め誤差が0.5ulp（丸める桁の0.5倍）にかなり近い状況は発生することがあります。
たとえばx = 0x1.000002p+0f（1.0f + std::numeric_limits<float>::epsilon()とかstd::nextafter(1.0f,2.0f)などとも表せます。1+2^-23です）な場合にstd::sqrt(x) == 1.0fになりますが、その丸め誤差は0.499999985ulpになります。
平方根を取ってから二乗するとどうなるか
平方根を計算してから、その値を二乗すると元に戻らないことがあります。
たとえば、x = 0x1.ffd508p+0fに対して、std::sqrt(x) == 0x1.69fab4p+0f（真の値は大体0x1.69fab4fffe986p+0くらい）であり、std::sqrt(x)*std::sqrt(x) == 0x1.ffd506p+0fになりますが、左辺の結果を丸める前は0x1.ffd5052c0e90p+0であり、真の値から1.413951ulpもずれた値になっています。
しかし、このずれはほぼ√2 ulpに収まります。つまり、1.5ulpを超えることはないので、xとstd::sqrt(x)*std::sqrt(x)の間の違いは最大でも1ulpに収まります。
二乗してから平方根を取るとどうなるか
なんと必ず元の値に戻ることが証明されているようです。
https://fixedpoint.jp/2016/03/18/sqrt-of-square-of-fp-number.html
もちろん、二乗した値が0になったり、INFINITYになった場合は元に戻ることはありません。
どういったアルゴリズムで計算できるのか
超越関数の場合、意地悪な引数を与えられると丸め誤差が0.5ulpに非常に近くなり、どちらに丸めてよいかを判断するために追加で計算する桁数が多くなる例があるようです。 平方根関数の場合、そのようなことはなく、仮数部の桁数から計算可能な定数桁を計算すれば丸める方向が決まるのでしょう（調べられませんでしたが、IEEEで平方根関数の丸めは正確に行うと決められているのですから、何かしらの裏付けがあるのだと思います）。
単純な方法としては、整数上で十分な精度で計算してそれを浮動小数点数に変換するというものがあります。 しかし、この方法は明らかに遅く、実用に耐えるものとは言えません。 なんとか、浮動小数点演算をうまく活用して平方根関数を実現できないものでしょうか。
以下のような実装をしてみたところ、全ての浮動小数点数の値でstd::sqrtと値が一致しました。ただし、std::sqrtと違い、エラー報告をしない点は動作が違いますので注意してください。
負の数を入れると-nanが返ってくるんですね。
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 invSqrt( 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;
}
float mySqrt( float x ) {
if( x == 0 ) {
// same sign
return x;
}
if( x < 0 ) {
return std::copysign( std::nanf(""), -1.0f );
}
if( std::isinf( x ) ) {
return x;
}
if( x < std::numeric_limits<float>::min() ) {
return std::ldexp( mySqrt( std::ldexp( x, 32 ) ), -16 );
}
const float q = invSqrt(x);
const float t = std::fma( x*q, x*q, -x );
const float X = x*q - std::ldexp( q*t, -1 );
return X;
}
apploxInvSqrt関数は、Fast Inverse Square Rootと呼ばれるビット演算テクニックを使っています（C++ではビット列を読み替えるためにキャストやunionを使うのはstrict alias rulesに反しますので、std::memcpyを使います。C++20ではstd::bit_castというそのものずばりの関数が導入されるのでそのような間違いは減ってゆくでしょう）。
invSqrt関数では、apploxInvSqrtで得られる初期値をもとに、ニュートン・ラフソン法を用いて逆数平方根の精度を上げています。もともと4bitは精度があるため、3回加速すれば8bit, 16bit, 32bitと精度が上がり、24bit精度のfloatでは十分な精度になります。
逆数平方根の精度が重要な場合は、最後の加速は丁寧に行う必要がありますが（参考：逆数平方根の収束公式と丸め誤差）、今回は平方根を求めることができればいいため、細かいことは気にしない実装になっています。
mySqrt関数では、基本的にはmySqrt(x) = x * invSqrt(x)の式によって計算を行います（これは、低速な除算を避けるためです。平方根自体をニュートン・ラフソン法で求めるのではなく、逆数平方根をニュートン・ラフソン法で求めたのも除算を避けるためです。平方根をニュートン・ラフソン法で求めたい場合、除算が必要になってしまいます）。
これだけだと最後の一桁が正しくならないことがあるので、（一次の）ニュートン法で補正します。微分係数は1/2√xであり、二次収束が必要な場合は1/2(x*q)の計算が必要になりますが、一次収束で十分なのでq/2で代用できます。
なお、非正規化数の平方根を求める場合、逆数平方根を介すと大きくなりすぎて表現できなくなる危険性があります。そのため、非正規化数であった場合は232をかけた値の平方根を求め、それを216で割ることにより、この問題を回避しています。
標準でない丸めモードを使うとどうなるのか
最近接丸め（丸め誤差がぴったり0.5ulpになることは平方根演算では発生しないので、偶数丸めでも無限大への丸めでもどちらでもよい）の場合は上記の方法で求められることが分かりました。
方向丸め（正の無限大への丸め、負の無限大への丸め、零への丸め）の場合もうまくいくのでしょうか？
結論から言うと、上記の方法では正しく求められません。
平方根関数の丸めモードを指定する場合、以下のようなコードになるでしょう。
std::fesetround(FE_TOWARDZERO); std::cout << std::sqrt(2.0f) << std::endl;
このように書いた場合、std::sqrt関数がライブラリ実装されており、その中で浮動小数点演算を行っている場合、その浮動小数点演算の丸めモードが意図せず変わってしまうことになります。
たとえば、零への丸めの場合、上記の方法ではx = 1.0fの時、q = invSqrt(x)についてq == 0x0.fffffep+0となってしまい、t = std::fma(x*q,x*q,-x)についてt == -1e-23になり、mySqrt(x) == 0x0.ffffffp+0となり、正しくありません。
tは正負両方の値を取りうるため、三つの丸めをすべて異なる方法で取り扱う必要が出てきます。いったいどうやったらよいのでしょう……？分岐してよいのなら、近似値を計算した後にstd::fmaで検算をして、ずれていれば1ulp増減させるという方法がありますが……。
やはりライブラリ実装をするのならば、整数上で計算するのが無難だということでしょうか。
しかしstd::sqrtは正しく計算するようです。しかも高速です。SSE命令のsqrtss命令を使っているからのようです。
そういえば先ほどの実装ではstd::fmaはちゃんと精度よく求まることを前提としていましたが、std::sqrtがライブラリ実装のアーキテクチャでstd::fmaだけはハードウェア命令で実現されるなんてことはあるのかと思われるかもしれません。実際、Intel Itanium (IA-64)では、FMA命令はあるものの、FDIVやFSQRT命令は存在せず、そのようなアーキテクチャの一例になっています。
