Rust で言うところの下記のような関数です。
fn is_square(nn: u64) -> bool { let n = (nn as f64).sqrt() as u64; n.wrapping_mul(n) == nn }
fn is_square(nn: u64) -> bool { ((nn as f64).sqrt() as u64).wrapping_pow(2) == nn }
これが問題ないということを示そうねというのが今回のお話です。
力技
u64 の範囲の平方数なんてのは $2^{32}$ 個しかないので、それらとその off-by-one を全部テストすることは数秒で可能です。
#[test] fn boundaries() { assert!(is_square(0)); assert!(is_square(1)); assert!(!is_square(2)); for i in 2..=u32::MAX as u64 { let ii = i * i; assert!(is_square(ii)); assert!(!is_square(ii - 1)); assert!(!is_square(ii + 1)); } assert!(!is_square(u64::MAX)); }
数学
数式での証明もしましょう。
IEEE 754 を前提として、sqrt は correctly-rounded で計算できるとします*1。
なので sqrt は問題ないのですが、整数型からの変換 (as f64) の部分で誤差が生じます。
それを考慮した上で、問題なく計算できるというのが面白いところですね。
double (binary64) を前提として、$p = 53$ とします。
Observation 1: $\roundp{\sqrt{94906267^2}} = 94906267$.
ここで、$94906265^2\lt 2^p\lt 94906267^2$ です。
>>> 94906267**2 9007199515875289 >>> 9007199515875289.0 9007199515875288.0 >>> math.sqrt(9007199515875289.0) 94906267.0
$$ \begin{aligned} \roundp{94906267^2} &= 94906267^2-1 \\ \sqrt{94906267^2-1} &= 94906266.{\footnotesize 999999994731644}{\scriptsize 012507624877103}{\tiny 348660069207956{\dots}} \\ \roundp{\sqrt{94906267^2-1}} &= 94906267. \end{aligned} $$
Lemma 2: 実数 $x\gt 0$ に対して、$1+\tfrac12x \gt \sqrt{1+x}$ が成り立つ。
Proof
$f(x) = (1+\tfrac12x)-(1+x)^{1/2}$ とする。$f(0) = 0$ なので、$x\gt 0 \implies f'(x)\gt 0$ を示せば十分。
両辺の差を微分して、 $$ \begin{aligned} f'(x) &= \tfrac12-\tfrac12(1+x)^{-1/2} \\ &= \tfrac12(1-(1+x)^{-1/2}). \end{aligned} $$ $x\gt 0$ のとき $(1+x)^{-1/2} \lt 1$ より従う。$\qed$
note (fact): $\sqrt{1+x} = 1 + {\small 0.5}x - {\small 0.125}x^2 + {\small 0.0625}x^3 - {\small 0.0390625}x^4 + O(x^5)$.
Lemma 3: 実数 $0\lt x\lt 1$ に対して、$1-x\lt\sqrt{1-x}$ が成り立つ。
Proof
$0\lt 1-x \lt 1$ であるから、$(1-x)^1 \lt (1-x)^{1/2}$ は明らか。$\qed$
Lemma 4: $f(x)$ を次で定義する。 $$ f(x) = \frac{\tfrac x2}{(1-\tfrac x2)-\sqrt{1-x}}. $$ このとき、任意の $0\lt x\lt 1$ に対して $f(x) \gt \tfrac 4x-3$ が成り立つ。
Proof
$$ \begin{aligned} f(x) &= \frac{\tfrac x2}{(1-\tfrac x2)-\sqrt{1-x}} \\ &= \frac{\tfrac x2}{(1-\tfrac x2)^2-(1-x)}\cdot( (1-\tfrac x2)+\sqrt{1-x}) \\ &= \frac{\tfrac x2}{\tfrac{x^2}4}\cdot( (1-\tfrac x2)+\sqrt{1-x}) \\ &\gt \tfrac2x\cdot( (1-\tfrac x2)+(1-x) ) \\ &= \tfrac2x\cdot(2-\tfrac32x) \\ &= \tfrac4x-3.\quad\qed \end{aligned} $$
note (fact): $f(x) = \tfrac4x-2-{\small 0.25}x-{\small 0.125}x^2-{\small 0.078125}x^3-{\small 0.0546875}x^4+O(x^5)$.
Therorem 5: 整数 $2^{p/2}\le n\le 2^p$ に対し、$\roundp{\sqrt{\roundp{n^2}}} = n$ が成り立つ。
Proof
この範囲において $\roundp{n}=n$ であることに注意する。 $\roundp{\sqrt{\roundp{n^2}}}\le n$ かつ $\roundp{\sqrt{\roundp{n^2}}}\ge n$ を示せばよい。
( $\le$ ): $\sqrt{\roundp{n^2}} \lt n + \hfloor{n}\cdot 2^{-p}$ を示す。
$$ \begin{aligned} \sqrt{\roundp{n^2}} &\le \sqrt{n^2(1+2^{-p})} \\ &= n\sqrt{1+2^{-p}} \\ &\lt n\cdot(1+2^{-(p+1)}) \\ &\lt n+\hfloor{n}\cdot 2^{-p}. \end{aligned} $$
( $\ge$ ):
Case 1: $\hfloor n=n \implies \sqrt{\roundp{n^2}} \ge n - \hfloor{n}\cdot 2^{-(p+1)}$ を示す。
$\hfloor n=n$ より、ある整数 $k$ に対して $n=2^k$ となる。 すなわち、$\roundp{n^2} = \roundp{2^{2k}} = 2^{2k}$ であり、$\sqrt{\roundp{n^2}} = n$ から、明らかに従う。
Case 2: $\hfloor n\lt n \implies \sqrt{\roundp{n^2}} \gt n - \hfloor{n}\cdot 2^{-p}$ を示す。
$$ \begin{aligned} \sqrt{\roundp{n^2}} &\ge \sqrt{n^2(1-2^{-p})} \\ &= n\sqrt{1-2^{-p}}. \end{aligned} $$
$n\le 2\hfloor{n}-1$ より $\hfloor{n}\ge \tfrac{n+1}2$ なので、 $$ \begin{aligned} n-\hfloor n\cdot2^{-p} &\le n-\tfrac{n+1}2\cdot 2^{-p} \\ &= n\cdot(1-2^{-(p+1)})-2^{-(p+1)}. \end{aligned} $$ よって、下記を示せば十分。 $$ n\sqrt{1-2^{-p}} - \left(n\cdot(1-2^{-(p+1)})-2^{-(p+1)}\right) \gt 0. $$
これは、 $$ n \le 2^p \lt 2^{p+2}-3 \lt f(2^{-p}) = \frac{2^{-(p+1)}}{(1-2^{-(p+1)})-\sqrt{1-2^{-p}}} $$ より従う。$\qed$
Corollary 6: 整数 $0\le n\le 2^p$ に対し、$\roundp{\sqrt{\roundp{n^2}}} = n$ が成り立つ。
Proof
$0\le n\lt 2^{p/2}$ のときは、$\roundp{n^2} = n^2$ より明らか。Theorem 5 より従う。
note: $n=2^p+1$ のときに成り立たないのは明らか。
>>> 2**53+1 9007199254740993 >>> sqrt((2**53+1)**2) 9007199254740992.0
wrapping に関して
ここ追記です。
n.wrapping_mul(n) の部分で変なことにならないことを示しておきます。
n.wrapping_mul(n) が wrap するのは n が$2^{32}$ 以上のときです。一方、入力 nn は $2^{64}-1$ 以下なので、n は $2^{32}$ 以下になります（$\sqrt{2^{64}-1}\le 2^{32}$ かつ $\roundp{2^{32}} = 2^{32}$ より）。
そのため、n が $2^{32}$ となる場合のみ考えればよいです。
n が $2^{32}$ のときは n.wrapping_mul(n) == 0 ですが、nn == 0 と両立することはないので、誤って true を返すことはありません。
また、nn が $2^{64}$ となることは（u64 の範囲から）ないので、誤って false を返すこともありません。$\qed$
あとがき
これで安心して sqrt を使って平方数判定をして生きていけます。
結果が所望の整数になってくれるため、「念のため round しておく」のような処理も不要であることがわかります。
？？「u64::isqrt は stabilize されてますよ笑」
ところで、Theorem 5 の ( $\ge$ ) の Case 2 を示すのに苦戦しました。以前自分で書いた記事を見て「なるほどな〜」となりました。衰えを感じます。悔しかったので少し違う方針で示しました。
fact の部分の Taylor 展開は sympy.series などで求めています。「こんな不等式成り立つの？」という気持ちになったときに、とりあえずグラフや Taylor 展開を見て「たしかに成り立ちそう・示せそうか〜」とやっていました。甘えかもしれません。
Gappa や Sollya や SageMath ともお友だちになろうとしているところです。
そのうちまた初等関数についての話を書くと思います。
おわり
おわりです。
