\(\sqrt{n}\) の計算に、浮動小数点型の sqrt を用いた際の誤差について記します。

\(\gdef\roundcirc#1{\operatorname{\circ}{(#1)}}\) \(\gdef\floor#1{\lfloor#1\rfloor}\) \(\gdef\ceil#1{\lceil#1\rceil}\)

IEEE 754 で定義される binary64 にあたる浮動小数点型で、丸めモード tiesToEven（表せる数のうち最も近いものに丸める。タイのときは仮数部が偶数になる方を選ぶ）で実数 \(x\) を丸めたときの値を \(\roundcirc{x}\) と書くことにします。

note: double (C++), float (Python), f64 (Rust) のデフォルトの挙動に相当します。

なお、C++ において sqrt(n) と書いたとき、n が long long だったとしても暗黙に n が double に変換されます。

例

下記の 2 数について考えます。

\[ \begin{aligned} n_1 &= (2^{26}+1)^2-1 \\ &= 67108865^2-1 \\ &= 4503599761588224, \\ n_2 &= \ceil{2^{26.5}}^2-1 \\ &= 94906266^2-1 \\ &= 9007199326062755. \end{aligned} \]

\(2^{52} \lt n_1 \lt 2^{53} \lt n_2\le 10^{18}\) が成り立ちます。

該当の型は仮数部のビット数が 53 なので、\(\roundcirc{n_1} = n_1\) が成り立ちます。また、\(\roundcirc{n_2} = n_2+1\) が成り立ちます。

\[ \sqrt{n_1} = 67108864.999999992549419{\dots} \]

であり、\(\roundcirc{\sqrt{n_1}} = \floor{\sqrt{n_1}}+1\) となります。 また、

\[ \sqrt{n_2+1} = 94906266 \]

が成り立ちます。こうして、\(\roundcirc{\sqrt{\roundcirc{n_1}}} = \floor{\sqrt{n_1}}+1\) および \(\roundcirc{\sqrt{\roundcirc{n_2}}} = \floor{\sqrt{n_2}}+1\) がわかりました。

\(n'=\roundcirc{n}\) の計算と、\(\roundcirc{\sqrt{n'}}\) の計算の 2 箇所で誤差が生じうるため、両方を気にする必要があります。\(n_1\) は後者、\(n_2\) は前者のみの誤差が生じる例です。

対策

誤差のことを気にしなくてよくなる簡単な方法としては、整数型で二分探索する方法が挙げられます。

また、sqrt(n) の返り値 m を前後 10 個くらいの i を調べ、i * i <= n && (i + 1) * (i + 1) > n となるものを探すという方法を使っている人もしばしばいます。実は、m と m - 1 についてのみ調べれば十分であることは示せます（解説ブログ）。

uint64_t uint64_floor_sqrt(uint64_t n) {
    if (n == 0) return 0; 
    uint64_t tmp_m1 = std::sqrt(n) - 1;
    return tmp_m1 * (tmp_m1 + 2) < n ? tmp_m1 + 1 : tmp_m1;
}

fn u64_floor_sqrt(n: u64) -> u64 {
    let tmp = (n as f64).sqrt() as u64;
    let tmp_m1 = tmp.saturating_sub(1);
    if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 }
}

また、Python では math.isqrt が用意されているので、それを使うのがよいでしょう。