平方数かどうか、改めdoubleの丸めに関する話

cpp

http://topcoder.g.hatena.ne.jp/nitoyon/20090208/1234068839

「c++の1.0の型って何だろう」って思った。まあsqrtはdoubleの引数を取れるみたいなのでdoubleで渡すということにしてみる。doubleの仮数部の精度は52bitだから、大雑把に言って12bit足りない。10進法で3桁半たりない。ということは10**18は入らないはず。

int main(){
  ULL n;
  cin >> n;
  DP(n * n);
  DP(ULL(double(n * n)));
  DP(sqrt(double(n * n)));
  ULL rt = ULL(sqrt(double(n * n)));
  DP(rt);
  DP(rt * rt);
}

4200000001
n * n: 17640000008400000001
ULL(double(n * n)): 17640000008400001024
sqrt(double(n * n)): 4.2e+09
rt: 4200000001
rt * rt: 17640000008400000001
4200000002
n * n: 17640000016800000004
ULL(double(n * n)): 17640000016800000000
sqrt(double(n * n)): 4.2e+09
rt: 4200000002
rt * rt: 17640000016800000004

あれ？2番目の例ではunsigned long longの17640000016800000004がdoubleの精度では表現できなくて少し小さな17640000016800000000になっているわけなので、sqrtを取った結果も真のsqrtである4200000002より少し小さな値になってunsigned long longに戻す際に正しい値よりも1小さい結果になるってのを期待したのにそうならないなぁ。僕は何を勘違いしているんだろう。


最初cinから取ってなかったんだけど、コンパイル時に計算されてしまってるかもと思ってcinから取るようにして、それでもこの挙動だから-O0をつけてコンパイルし直してみたりしたんだけどやっぱり変わらないなぁ。

>>> 4200000001 * 4200000001
17640000008400000001
>>> from math import log
>>> log(_) / log(2)
63.935484364442317

うーん、うっかり一桁足りないとかじゃないなぁ。ちゃんと64bitギリギリな値になっている。4200000001から上に順に探していったけど正しくない値を返す最初の数は

n: 4294967296
n * n: 0

だなぁ。


なるほど。FLT_ROUNDSの値を表示してみたら「(__builtin_flt_rounds ()): 1」と出力された。1は「最も近い数へ丸める」を意味する。(ref. Manpage of FENV)

ISO C では、丸めモードを最も近い値への丸め (FE_TONEAREST) に設定し、すべての例外をクリアし、不停止 (non-stop) (例外が起きても継続する) モードとするように規定されている。

つまり僕の「doubleをunsigned long longに変換する際の丸めは切り捨てである」という理解が間違っていたということだな。

# 「最も近い整数」が複数あるケースの挙動が言及されていないように見えるけどまぁ本題とずれるから今日はこれくらいにしとこう
修正：違う違う、全然違った。ここで書かれているのはdoubleを整数に丸める話ではなくてsqrtの正確な計算結果を正確な値を扱えないとびとびのdoubleで表現する際にどう丸めるかという話ね。

というわけで本題に戻ると、「その方法で平方数かどうか判定すると52bit以上の時に微妙に小さくなって切り捨てで下の整数になってしまって正しく判定できない！」というツッコミは空振りでした。

しかし結局「doubleに変換してsqrtして2乗したものが元の値かどうかで平方数かどうか判定」をしていいのか、やっぱりしてはいけないのかがわからない……