符号なしの 64-bit 整数型の値に関して平方根を計算したい局面はそこそこあります。 しかし、多くの言語ではそのための関数が用意されておらず、浮動小数点数型用のメソッドを使う人が多くいます。
今回は、そうした方針が妥当か？（誤差が出うることを踏まえつつ、正しい値を復元するのは簡単か？）という話です。
記法・前提
浮動小数点数型に対する一般的な知識はあるとします。ない人向けの記事はそのうち書きます。
実数 $x$ を該当の浮動小数点数型の値に丸めたものを $\roundp x$ と書くことにします。
IEEE 754 の規格において、sqrt(x) は $\roundp{\sqrt x\,}$ を返すことを規定しているので、ここではそれを仮定します。
note: 他の関数の誤差に関して
GCC において、sin(x), log(x), cbrt(x) などは、$\roundp{\sin(x)}$, $\roundp{\ln(x)}$, $\roundp{\sqrt[3]{x}\,}$ を返すとは限りません*1。sqrt などは IEEE 754 の規格で required の枠に含まれるのに対し、sin などは recommended の枠に含まれるので、建前上は「recommended の方は提供していませんよ〜」ということなのでしょうか。今後、正確に丸めた (correctly-rounded) $\roundp{\sin(x)}$ を返す crsin(x) のような関数を提供するかも？とのことです。
任意の正規化数 $x$ に対し、$2^e \le x\lt 2^{e+1}$ が成り立つような $e$ が一意に存在します。この $e$ を用いて、$\hfloor x=2^e$ と定義します。 $\hfloor x\le x\lt 2\hfloor x$ や $\tfrac x2\lt\hfloor x\le x$ が成り立ちます。 また、$x$ が正整数の下では $\hfloor x\le x\le 2\hfloor x-1$ や $\tfrac{x+1}2\le\hfloor x\le x$ が成り立ちます。
下記では、次のような関数を作ることを想定します。
/// $\floor{\sqrt n}$ を返す。 fn uint_floor_sqrt(u64 n) -> u64 { f64 tmp = (n as f64).sqrt(); // tmp が持ちうる誤差に対して何かしらの補正をする todo!() }
どのように補正を行うかについてはまだ与えないことにします（適当に考えた方針が「実はよかったよね」とするのではなく、考察ベースで自然なアルゴリズムを考えたいため）。考察次第では、.sqrt() をする前に何かしらをするかもしれません。
？？「仮数部が 64 bits 以上の
long doubleがある環境ならlong doubleでやればよくないか？」
それはそう... そうした型のない言語を想定しているとかの建前でやり過ごします。
考察
note: $n \le 2^{53}$ では $\roundp n=n$ ですが、それより大きい場合はそうとは限らないため、気をつける必要があります。
$\gdef\roundsqrt#1{\roundp{\sqrt{#1}\,}}$
Observation 1:
u32 に関してはすぐに全探索できるため、その範囲を f32（仮数部 24 bits）で計算してみます。
fn main() { for i in 0..=u16::MAX as u32 { for j in i * i..=i * i + 2 * i { let r = (j as f32).sqrt() as u32; assert_eq!(r, i, "failed on {}", j); } } println!("passed all tests."); }
assertion `left == right` failed: failed on 16785407 left: 4097 right: 4096
やはり、これはどうやらうまくいかないようです。 $$ \floor{\roundsqrt{(2^{12}+1)^2-2}} = 2^{12}+1. $$
note: f64 においても、$\floor{\sqrt n}$ を固定して二分探索を使ったり、固定した際の最大値を試してから降順に見ていくなどで、同様の探索ができます。
最小の反例は 4503599761588224 で、下記のようになっています。
$$
\floor{\roundsqrt{(2^{26}+1)^2-1}} = 2^{26}+1.
$$
$(2^{\floor{p/2}}+1)^2$ より少し小さいあたりでだめになってしまう感じでしょうか。
Observation 2:
では、$+1$ の誤差は許容してみましょう。
fn main() { for i in 0..=u16::MAX as u32 { for j in i * i..=i * i + 2 * i { let r = (j as f32).sqrt() as u32; assert!([i, i + 1].contains(&r), "failed on {}", j); } } println!("passed all tests."); }
passed all tests.
これはよいらしいです。
Observation 3:
ということで、u64 にして次のような関数を考えてみます。
fn u64_floor_sqrt(n: u64) -> u64 { let tmp = (n as f64).sqrt() as u64; let tmp_m1 = tmp.saturating.sub(1); // if tmp >= 1 { tmp - 1 } else { tmp }; // ((tmp - 1) + 1)**2 <= n, in an overflow-safe way if tmp_m1 * (tmp_m1 + 2) < n { tmp } else { tmp_m1 } }
$n = 0$ の場合や、$n = 2^{64}-1$ で $2^{32}$ が返ってきた場合などにオーバーフローしないように気をつけています。 そうしたケースに気を払う必要がない状況では、もう少し雑にやっても大丈夫そうです。
fn main() { for i in 0..=u32::MAX as u64 { let actual_left = u64_floor_sqrt(i * i); let actual_right = u64_floor_sqrt(i * i + 2 * i); let expected = i; assert_eq!(actual_left, expected); assert_eq!(actual_right, expected); } println!("passed all tests."); }
passed all tests.
これもよいらしいです。$\floor{\sqrt n}-1$ や $\floor{\sqrt n}+2$ になることはなさそうです。
$\roundp x$ や $\sqrt x$ の単調性から、u64_floor_sqrt は正当そうです。
これで「u64 の範囲ではこれでいけます、めでたし」で終わってもいいにはいいんですが、たとえば C++ において unsigned __int128 の範囲を __float128 で計算するのを同じ方針で調べることはできないため、もう少し数学的に考えます。
もしかすると 128-bit では反例が見つかるかもしれないし、そうではないかもしれませんから。
ここまでの話をまとめると、次のようなことを示したいです。
Conjecture 4: それなりに大きい $n$ に対しても $\floor{\sqrt n} \le \roundsqrt{\roundp n} \lt \floor{\sqrt n}+2$ が成り立つ。
まず、浮動小数点数 $\hat x$ と実数 $x$ に対して、$\roundp x = \hat x$ となる十分条件について触れておきます。 $$ \roundp x=\hat x \impliedby \begin{cases} \hat x-\hat x\cdot 2^{-1-p} \le x\le \hat x+\hat x\cdot 2^{-p}, & \text{if }\hfloor{\hat x}=\hat x; \\ \hat x-\hfloor{\hat x}\cdot 2^{-p}\lt x\lt\hat x+\hfloor{\hat x}\cdot 2^{-p}, & \text{if }\hfloor{\hat x}\gt\hat x. \end{cases} $$ $\roundp x\ge \hat x$ などを示したいときは片側の不等式が成り立てばよいです。
$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、整数 $n\in[k^2\ldots(k+1)^2)$ が conjecture 4 を満たすことを示します。
整数性から $n\in[k^2\ldots k(k+2)]$ であることと、単調性から $$\roundsqrt{\roundp{k^2}}\le\roundsqrt{\roundp n}\le \roundsqrt{\roundp{k(k+2)}}$$ であることに注意します。そのため、 $$ k\le \roundsqrt{\roundp{k^2}} \wedge k+2\gt \roundsqrt{\roundp{k(k+2)}} $$ を示せば十分です。考えている状況に合わせ、$\roundp{k^2}\lt\infty$ である前提とします。 また、この範囲で $\roundp k=k$ が成り立つことに注意しておきます。 加えて、$p\ge 3$ も仮定しておくことにします。
Lemma 5: $0\le k\le 2^p$ なる整数 $k$ に対して、$k\le \roundsqrt{\roundp{k^2}}$ が成り立つ。
disclaimer: もともと $k\le 2^p$ で示す予定だったのですが、Lemma 6 が $k\le 2^{p-1}$ でしか成り立たなさそうということが Lemma 5 を示した後にわかりました。Lemma 5 を直すのが面倒だったのでそのままになっています。
Proof
ある $e$ が存在して $k = 2^e$ のとき、$2^{2e} = 2^{p-1}\times 2^{2e+1-p}$ と表せるため、$\roundp{2^{2e}}=2^{2e}$ となる。 $\roundsqrt{\roundp{k^2}} = \roundsqrt{k^2} = \roundp k=k$ が成り立つ。すなわち、$k\le \roundsqrt{\roundp{k^2}}$ となる。
$k = 0$ のときも同様にして成り立つことが示せる。
次に、$k = 2^e$ なる $e$ が存在しないかつ $k\ge1$ とし、 $$ k-\hfloor{k}\cdot 2^{-p}\lt \sqrt{\roundp{k^2}} $$ を示す。これらを評価すると、 $$ \begin{aligned} k-\hfloor k\cdot 2^{-p} &\le k-\tfrac{k+1}2\cdot 2^{-p} \\ &= (1-2^{-1-p})k - 2^{-1-p}; \\ \sqrt{\roundp{k^2}} &\ge \sqrt{k^2-\hfloor{k^2}\cdot 2^{-p}} \\ &\ge \sqrt{k^2-k^2\cdot 2^{-p}} \\ &= k\cdot \sqrt{1-2^{-p}} \end{aligned} $$ となるので、$(1-2^{-1-p})k-2^{-1-p}\lt k\cdot\sqrt{1-2^{-p}}$ を示せばよい。非負であることは明らかであるから二乗の差を考え、 $$ \begin{aligned} &\phantom{{}={}} (k\cdot\sqrt{1-2^{-p}})^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ &= (1-2^{-p})k^2 - ( (1-2^{-1-p})k-2^{-1-p})^2 \\ %&= (1-2^{-p})k^2 - ( (1-2^{-1-p})^2k^2-2^{-p}\cdot(1-2^{-1-p})k+2^{-2-2p}) \\ %&= (1-2^{-p})k^2 - ( (1-2^{-p}+2^{-2-2p})k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -(2^{-2-2p}\cdot k^2-2^{-p}\cdot (1-2^{-1-p})k+2^{-2-2p}) \\ %&= -2^{-2-2p}\cdot k^2+2^{-p}\cdot (1-2^{-1-p})k-2^{-2-2p} \\ %&= 2^{-p}\cdot (-2^{-2-p}\cdot k^2 + (1-2^{-1-p})k - 2^{-2-p}) \\ %&= 2^{-2p}\cdot (-2^{-2}\cdot k^2 + (2^p-2^{-1})k - 2^{-2}) \\ &= 2^{-2-2p}\cdot (-k^2 + (2^{2+p}-2)k - 1) \\ \end{aligned} $$ となる。
$-k^2+(2^{2+p}-2)k-1$ の部分が、主要項が負の $2$ 次関数になっていることに注意すると、これが $k=1$ および $k=2^p$ で正になっていることを示せば、$1\le k\le 2^p$ の範囲で正であることがわかる。 $$ \begin{aligned} -1^2 + (2^{2+p}-2)\cdot 1-1 &= 4\cdot(2^p-1); \\ -(2^p)^2 + (2^{2+p}-2)\cdot 2^p-1 %&= -2^{2p} + 2^{2+2p} - 2^{1+p} - 1 &= 2^p\cdot(3\cdot 2^p-2)-1 \end{aligned} $$ より、$p\ge 1$ でこれらは正となることがわかる。
以上より、$0\le k\le 2^p$ なる任意の整数 $k$ に対して $k\le\roundsqrt{\roundp{k^2}}$ が成り立つ。$\qed$
Lemma 6: $0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。
Proof
任意の $0\le k\le 2^{p-1}$ に対して、 $$ \sqrt{\roundp{k(k+2)}} \lt (k+2)-\hfloor{k+2}\cdot 2^{-p} $$ を示せばよい。両辺を評価して $$ \begin{aligned} \sqrt{\roundp{k(k+2)}} &\le \sqrt{k(k+2)+\hfloor{k(k+2)}\cdot 2^{-p}} \\ &\le \sqrt{k(k+2)+k(k+2)\cdot 2^{-p}} \\ &= \sqrt{k(k+2)}\cdot \sqrt{1+2^{-p}}; \\ (k+2)-\hfloor{k+2}\cdot 2^{-p} &\ge (k+2)-(k+2)\cdot 2^{-p} \\ &= (k+2)\cdot (1-2^{-p}) \end{aligned} $$ となるので、 $$ \sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}} \lt (k+2)\cdot(1-2^{-p}) $$ を示す。Lemma 5 同様にして、 $$ \begin{aligned} &\phantom{{}={}} ( (k+2)\cdot(1-2^{-p}))^2 - (\sqrt{k(k+2)}\cdot\sqrt{1+2^{-p}})^2 \\ &= (k+2)^2\cdot(1-2^{1-p}+2^{-2p}) - k(k+2)\cdot(1+2^{-p}) \\ &= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)(k+2)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (2^{2p}-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-(2^{2p}+2^p)k\right) \\ %&= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}+1)k+2(2^{2p}-2^{1+p}+1)-2^p\cdot k\right) \\ &= 2^{-2p}(k+2)\cdot\left( (-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)\right) \\ \end{aligned} $$ となる。$(-2^{1+p}-2^p+1)k+2(2^{2p}-2^{1+p}+1)$ の部分が、主要項が負の $1$ 次関数であることに注意すると、これが $k=0$ および $k=2^{p-1}$ で正であることを示せば、$0\le k\le 2^{p-1}$ の範囲でこれが正であることがわかる。$p\ge 3$ の前提で、 $$ \begin{aligned} &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 0+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}-2^{2+p}+1; \\ &\phantom{{}={}} (-2^{1+p}-2^p+1)\cdot 2^{p-1}+2(2^{2p}-2^{1+p}+1) \\ &= 2^{1+2p}+2^{p-1} - (2^{2p}+2^{2p-1}+2^{2+p}) + 2 \\ &\ge 2^{1+2p}+2^{p-1} - 2^{1+2p} + 2 \\ &\gt 2^{p-1}+2 \end{aligned} $$ より、これらは正であることがわかる。
すなわち、$0\le k\le 2^{p-1}$ なる整数 $k$ に対して、$k+2\gt \roundsqrt{\roundp{k(k+2)}}$ が成り立つ。$\qed$
よって、Lemmata 5–6 より、$0\le n\le 2^{2(p-1)}$ なる整数 $n$ に対して $$\floor{\sqrt{n}} \le \floor{\roundsqrt{\roundp n}} \le \floor{\sqrt{n}}+1$$ が成り立ちます。すなわち、上記の関数の正当性が示されたことになります。
note: (n as f64).sqrt() as u64 で計算しているのが $\floor{\roundsqrt{\roundp n}}$ です。n as f64 が $\roundp n$ に、x.sqrt() が $\roundsqrt x$ に、r as u64 が $\floor r$ に対応します。
実装例
Rust
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 } }
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; }
def int_floor_sqrt(n: int) -> int: tmp = int(n**0.5) return tmp if tmp**2 <= n else tmp - 1
下記のケースでテストしています。
| $n$ | $\floor{\sqrt n}$ |
|---|---|
| $0$ | $0$ |
| $1$ | $1$ |
| $(2^{26}+1)^2-1$ | $2^{26}$ |
| $(2^{26}+1)^2$ | $2^{26}+1$ |
| $2^{64}-1$ | $2^{32}-1$ |
あとがき
つかれました。おなかがすきました。
証明は一旦 iPad で下書きをしているのですが、組版する段階になって符号が間違っていることに気づいたり、指数部の符号も間違っていたり、示したい命題が違っていたりで散々でした。 下書きの証明をそのまま写すのではなく、組版しながらまた別途考えているので、「あれ、こんな式書いた記憶なくない？」というところで気づいたりしています。
さておき、ある程度は浮動小数点数とも仲よくなってきたんじゃないかなという気持ちがあります。
初心者の頃、「平方根を sqrt(n) で求めるのは誤差が出ちゃう」と聞いて、返り値の ±30 くらい前後を探すようなコードを書いていたのが懐かしいです。
また、今回の話を踏まえると、素数判定のときにループの上限を i <= sqrt(n) とするのは正当そうということもわかりますね。
約数列挙についてはどうでしょう。あまり考えていませんが、少なくとも u64・f64 で全探索した限りでは、$n=i\cdot(i+1)$ のとき $\floor{\roundsqrt{\roundp n}}=i$ となるようだったので、重複して列挙されることはなさそうな気がします。
いかがでしたでしょうか。
今回の話は、ここ数回の中でも競プロ文脈に比較的近そうな気がします。
とはいえ、この記事を見て「なるほど！」となって、rated のコンテストでも f64 ベースの実装で平方根を求める人がいたらなかなか勇気があるなという気もします。
そんなこともないかも。
証明が嘘だったときに責任を取るのは嫌なので、各自で証明し直してほしい気はします。
当初はこの記事は「書く予定リスト」に含まれていなかったため、予告していた「浮動小数点数初心者向け」の記事は次の次になりそうです。ごめんね。
おわり
おわりです。
*1:「有限桁なので、数学的な厳密な値からは誤差が出て当然」という話ではなく、「無限精度で計算した値を正確に丸めた値が欲しいが、それとは異なる値が返ってきうる」という話です。
