ABC で等比数列の判定をしたくなることはしばしばあります。
これに対して double で判定するとうまくいかないケースがあり、たとえば下記です。
$$
67114657 \oslash 67114656 = 67114658 \oslash 67114657 = \texttt{1}.\texttt{0000003FFE96}_{(16)}.
$$
すなわち、$\angled{67114656, 67114657, 67114658}$ を等比数列であると判定してしまいます。
今回は、仮数部が $53$-bit ではなく $64$-bit や $113$-bit であれば正当なのか？というのを考えていきたいです。
記法
$% \gdef\poslash{\stackrel{p}{\oslash}}$ $\gdef\poslash{\oslash_p}$
$A_k = \{n\in\Z\mid 0\lt n\le k\}$ とします。
仮数部が $p$-bit の場合の除算を $\poslash$ と書くことにします。丸めモードは tiesToEven です。
$p=53$, $p=64$, $p=113$ の型を、それぞれ f64, f80, f128 と呼ぶことにします。AtCoder 環境においては、それぞれ double, long double, __float128 に相当します。
前提
ただし、数列 $S = (S_1, S_2, \dots, S_N)$ が等比数列であるとは、ある実数 $r$ が存在してすべての整数 $1\le i\lt N$ に対して $S_{i+1} = rS_i$ が成り立つことをいいます。
この定義に沿うと、たとえば $(1, 0)$ は等比数列であり、$(0, 1)$ は等比数列ではないことに注意しましょう。以下では $S_i \ne 0$ を仮定します*1。
これに沿った素朴な実装は次のようなものだと思います。
template <typename Float> bool is_geometric_nonzero(const std::vector<Float>& s) { if (s.size() <= 1) return true; float_type r = s[1] / s[0]; for (size_t i = 0; i < n - 1; ++i) { if (s[i + 1] != r * s[i]) return false } return true; }
これは、std::vector<f80>{7, 15} や std::vector<f128>{7, 29} に対して false を返します。よって、正当ではありません。
長さが $2$ のときは自明だとして if (s.size() <= 2) return true; に変えたとしても、std::vector<f80>{25, 15, 9} や std::vector<f128>{25, 55, 121} などが反例になります。
誤差を含みうる値をその後の計算に使うのを避けた次のコードを考えます。
template <typename Float> bool is_geometric_nonzero(const std::vector<Float>& s) { if (s.size() <= 1) return true; float_type r = s[1] / s[0]; for (size_t i = 0; i < n - 1; ++i) { if (s[i + 1] / s[i] != r) return false } return true; }
これに基づく解法が正当なのか？ということを調べていきます。 なお、ある実数 $r$ に対して $S_2 = rS_1 = r^2S_0$ が成り立つときは $S_2\oslash S_1 = S_1\oslash S_0 = \roundcirc{r}$ となるため、Yes のところを誤って No とすることはありません*2。以下では、No のところを誤って Yes とする場合について調べます。
本題
$x\oslash y = y\oslash z$ を満たすような組 $(x, y, z)\in A_{10^9}^3$ が存在するか？というのを考えていきます。 想定解は $xz = y^2$ を整数演算で行うものでしょうから、AtCoder あるある制約として $10^9$ を上限とするのは妥当な仮定だと思います。
Lemma 1: $k\ge 2$ のとき、任意の $y_1, y_2 \in A_k$ に対して $\lcm(y_1, y_2)\le k(k-1)$ が成り立つ。
Proof
Case 1: $(y_1, y_2) = (k, k)$
このとき $\lcm(y_1, y_2) = k \le k(k-1)$ が成り立つ。
Case 2: $(y_1, y_2) \ne (k, k)$
このとき $\lcm(y_1, y_2) \le y_1y_2 \le k(k-1)$ が成り立つ。$\qed$
Claim 2: $k\ge 2$ のとき、任意の $(x_1, y_1, x_2, y_2)\in A_k^4$ に対して、$x_1 / y_1 \ne x_2 / y_2$ ならば下記が成り立つ。 $$ \left| \frac{x_1}{y_1} - \frac{x_2}{y_2} \right| \ge \frac1{k(k-1)}. $$
Proof
一般性を失わず $x_1 / y_1 \gt x_2 / y_2$ とし、 $$ \frac{x_1}{y_1} - \frac{x_2}{y_2} \ge \frac1{k(k-1)} $$ を示す。
$g = \gcd(y_1, y_2)$ とし、$(y_1', y_2') = (y_1/g, y_2/g)$ とする。このとき、 $$ \begin{aligned} \frac{x_1}{y_1} - \frac{x_2}{y_2} &= \frac{x_1y_2 - x_2y_1}{y_1y_2} \\ &= \frac{x_1y_2' - x_2y_1'}{y_1'y_2' g} \\ &\ge \frac1{\lcm(y_1, y_2)} \\ &\ge \frac1{k(k-1)}. \quad\qed \end{aligned} $$
$x\oslash y \stackrel{?}{=} y\oslash z$ ではなく、より一般に $x_1\oslash y_1 \stackrel{?}{=} x_2\oslash y_2$ について考えます。
Claim 3: $k\le 2^{(p-1)/3}$ とする。任意の $(x_1, y_1, x_2, y_2)\in A_k^4$ に対して $x_1/y_1 \ne x_2/y_2$ ならば $x_1\oslash_p y_1 \ne x_2\oslash_p y_2$ が成り立つ。
Proof
ある $|\varepsilon_1|, |\varepsilon_2| \le 2^{-p}$ に対して $$ \begin{aligned} x_1\oslash_p y_1 &= x_1/y_1\cdot (1+\varepsilon_1), \\ x_2\oslash_p y_2 &= x_2/y_2\cdot (1+\varepsilon_2) \end{aligned} $$ が成り立つ。
一般性を失わず $x_1 / y_1 \gt x_2 / y_2$ とする。このとき、$x_1\oslash_p y_1 \ge x_2\oslash_p y_2$ は明らか。 $$ \begin{aligned} (x_1\oslash_p y_1)-(x_2\oslash_p y_2) &= x_1/y_1\cdot (1+\varepsilon_1)-x_2/y_2\cdot (1+\varepsilon_2) \\ &\ge x_1/y_1\cdot (1-2^{-p})-x_2/y_2\cdot (1+2^{-p}) \\ &= (x_1 / y_1 - x_2 / y_2) - 2^{-p}\cdot (x_1 / y_1 + x_2 / y_2) \\ &\ge \frac1{k(k-1)} - 2^{-p}\cdot 2k \\ &\gt k^{-2} - 2^{-(p-1)}\cdot k \\ &\ge 0 \end{aligned} $$ より $x_1\oslash_p y_1 \gt x_2\oslash y_2$ が従う。$\qed$
Corollary 4: 任意の $(x_1, y_1, x_2, y_2)\in A_{10^9}^4$ に対して $x_1/y_1 \ne x_2/y_2$ ならば $x_1\oslash_{113} y_1 \ne x_2\oslash_{113} y_2$ が成り立つ。
Proof. $10^9\lt 2^{(113-1)/3}$ より従う。$\qed$
すなわち、$|A_i|\le 10^9$ に対しては __float128 を用いた解法は正当ということですね。long double の解法についても考えていきます。
Claim 3 の証明において、 $$(x_1 / y_1 - x_2 / y_2) - 2^{-p}\cdot (x_1 / y_1 + x_2 / y_2) \ge \frac1{k(k-1)} - 2^{-p}\cdot 2k $$ としているのは大雑把すぎます。左側の項は $y_1, y_2\approx k$ としているのに右側の項は $y_1, y_2\approx 1$ としていて、全体としてはもっと厳しく評価できそうです。
整数 $p\ge 2$ に対し、$i\oslash_{p} (i-1) = (i+1)\oslash_{p} i$ なる最小の $i$ を $i_p$ とします。
Observation 5: $i_{64} = 3037041293$ であり、 $$ \begin{aligned} 3037041293 \oslash_{64} 3037041292 &= 3037041294 \oslash_{64} 3037041293 \\ &= \texttt{1}.\texttt{000000016A08A7B8}_{(16)} \end{aligned} $$ が成り立つ。$\eod$
Observation 6: $i_{53} = 67114657$ であり、 $$ \left(\frac{i_{64}}{i_{53}}\right)^2 = \left(\frac{3037041293}{67114657}\right)^2 \approx 2047.701 \approx 2^{64-53} $$ が成り立つ。$\eod$
Observation 7: $p$ が $1$ 増えるごとに $i_p$ はだいたい $\sqrt 2$ 倍になっていて、$i_p\approx 2^{(p-1)/2} + 2^{(p-3)/4}$ 程度である。$\eod$
全探索に基づく表
| $p$ | $i_p$ |
|---|---|
| $2$ | $3$ |
| $3$ | $4$ |
| $4$ | $5$ |
| $5$ | $6$ |
| $6$ | $7$ |
| $7$ | $11$ |
| $8$ | $15$ |
| $9$ | $20$ |
| $10$ | $26$ |
| $11$ | $37$ |
| $12$ | $53$ |
| $13$ | $70$ |
| $14$ | $98$ |
| $15$ | $137$ |
| $16$ | $191$ |
| $17$ | $268$ |
| $18$ | $375$ |
| $19$ | $529$ |
| $20$ | $741$ |
| $21$ | $1047$ |
| $22$ | $1465$ |
| $23$ | $2081$ |
| $24$ | $2948$ |
| $25$ | $4142$ |
| $26$ | $5832$ |
| $27$ | $8257$ |
| $28$ | $11603$ |
| $29$ | $16475$ |
| $30$ | $23284$ |
| $31$ | $32897$ |
| $32$ | $46508$ |
| $33$ | $65718$ |
| $34$ | $92937$ |
| $35$ | $131329$ |
| $36$ | $185773$ |
| $37$ | $262507$ |
| $38$ | $371062$ |
| $39$ | $524801$ |
| $40$ | $741728$ |
| $41$ | $1049301$ |
| $42$ | $1483930$ |
| $43$ | $2098177$ |
| $44$ | $2967454$ |
| $45$ | $4195753$ |
| $46$ | $5932969$ |
| $47$ | $8390657$ |
| $48$ | $11864339$ |
| $49$ | $16780113$ |
| $50$ | $23730607$ |
| $51$ | $33558529$ |
| $52$ | $47459580$ |
| $53$ | $67114657$ |
| $54$ | $94911152$ |
| $55$ | $134225921$ |
| $56$ | $189813286$ |
| $57$ | $268447042$ |
| $58$ | $379638923$ |
| $59$ | $536887297$ |
| $60$ | $759269842$ |
| $61$ | $1073764995$ |
| $62$ | $1518528457$ |
| $63$ | $2147516417$ |
| $64$ | $3037041293$ |
Observation 8: 下記が成り立つ。 $$ \begin{aligned} &\phantom{{}={}} 72057594227740468 \oslash_{113} 72057594227740467 \\ &= 72057594227740469 \oslash_{113} 72057594227740468 \\ &= \texttt{1}.\texttt{00000000000000FFFFFFF4AFB0CD}_{(16)}. \quad\eod \end{aligned} $$
「任意の $(x_1, y_1, x_2, y_2)\in A_k^4$ に対して、$x_1/y_1\gt x_2/y_2$ ならば $(x_1 / y_1 - x_2 / y_2) - r_k\cdot (x_1 / y_1 + x_2 / y_2) \ge 0$」を満たすような $r_k$ を考えたいです。
$$ r_k\le \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} $$ なので、この右辺を最小化すればよく、小さい $k$ について全探索で求めると下記のようになりました。 $$ (1/r_k)_{k=3}^{16} = (7, 17, 31, 49, 71, 97, 127, 161, 199, 241, 287, 337, 391, 449). $$
どうやら、$(x_1, y_1, x_2, y_2) = (k-1, k-2, k, k-1)$ で達成できそうです。 $$ \begin{aligned} \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} &\ge \frac{(k-1)(k-1)- k(k-2)}{(k-1)(k-1) + k(k-2)} \\ % &= \frac{(k^2-2k+1)-(k^2-2k)}{(k^2-2k+1)+(k^2-2k)} \\ &= \frac{1}{2k^2-4k+1} \\ \end{aligned} $$ が言えるのであれば、$r_k \le (2k_2-4k+1)^{-1}$ が言えるので、Claim 3 より厳しい評価ができますね。
Claim 9: $k\ge 3$ とする。任意の $(x_1, y_1, x_2, y_2)\in A_k^4$ に対して、$x_1/y_1\gt x_2/y_2$ ならば $$ \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} \ge \frac{1}{2k^2-4k+1} $$ が成り立つ。
Proof
$$ \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} \lt \frac{1}{2k^2-4k+1} $$ なる組 $(x_1, y_1, x_2, y_2)$ が存在しないことを示す。
Case 1: $x_1 = y_1$
$$ \begin{aligned} \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} &= \frac{y_2 - y_1}{y_2 + y_1} \\ &= \frac{1 - y_1/y_2}{1 + y_1/y_2} \\ \end{aligned} $$ となる。$y_2\gt y_1\gt 0$ より $0\lt y_1/y_2\lt 1$ であり、$f(x) = \frac{1-x}{1+x}$ が $0\lt x\lt 1$ で単調減少であることから、$\frac{y_1}{y_2} = \frac{k-1}k$ のとき最小となる。 $$ \begin{aligned} \frac{y_2 - y_1}{y_2 + y_1} &\ge \frac{k - (k-1)}{k + (k-1)} \\ &= \frac{1}{2k-1} \\ &\gt \frac{1}{2k^2-4k+1}. \end{aligned} $$
Case 2: $(x_2 = y_2) \vee (x_1 = x_2) \vee (y_1 = y_2)$
対称性より Case 1 と同様にして示せる。
Case 3: $(x_1, y_2) = (k, k)$
Cases 1–2 より、$(x_1\ne y_1) \wedge (x_2 \ne y_2) \wedge (x_1 \ne x_2) \wedge (y_1 \ne y_2)$ としてよく、$x_2y_1\le (k-1)^2$ が成り立つ。 Case 1 同様に $x_2y_1 = (k-1)^2$ のときに最小となるから、 $$ \begin{aligned} \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} &\ge \frac{k^2-(k-1)^2}{k^2 + (k-1)^2} \\ % &= \frac{k^2-(k^2-2k+1)}{k^2 + (k^2-2k+1)} \\ &= \frac{2k-1}{2k^2-2k+1} \\ &\gt \frac{1}{2k^2-4k+1}. \end{aligned} $$
Case 4: $\{x_1, y_2\} = \{k, k-1\}$
Case 3 同様にして、$x_2y_1\le (k-2)(k-3)$ が成り立つ。$x_2y_1 = (k-2)(k-3)$ のときに最小となるから、 $$ \begin{aligned} \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} &\ge \frac{k(k-1)-(k-2)(k-3)}{k(k-1)+(k-2)(k-3)} \\ &= \frac{4k-6}{2k^2-6k+6} \\ &\gt \frac{1}{2k^2-4k+1}. \end{aligned} $$
Case 5: ある $l\le k-2$ に対して $\{x_1, y_2\} = \{k, k-2\}$
Case 3 同様にして $x_2y_1\le (k-1)(k-2)$ が成り立つため、 $$ \begin{aligned} x_1y_2 + x_2y_1 &\le k(k-2) + (k-1)(k-2) \\ &= 2k^2-5k+2 \\ &\lt 2k^2-4k+1 \end{aligned} $$ である。また、$x_1y_2-x_2y_1\ge 1$ であるから、 $$ \begin{aligned} \frac{x_1y_2 - x_2y_1}{x_1y_2 + x_2y_1} &\gt \frac{1}{2k^2-4k+1}. \end{aligned} $$
Case 6: ある $l_1\le k-1$ および $l_2\le k-1$ に対して $\{x_1, y_2\} = \{l_1, l_2\}$
$x_1y_2 \gt x_2y_1$ より $(k-1)^2 \gt x_2y_1$ が成り立つ。よって、 $$ \begin{aligned} x_1y_2 + x_2y_1 &\le (k-1)^2 + (k-1)^2 - 1 \\ &= 2k^2-4k+1 \end{aligned} $$ であるから、Case 5 同様にして従う。$\qed$
よって、Claim 3 の証明の方針と合わせて、 $$ 2^{-p} \le \frac{1}{2k^2-4k+1} $$ すなわち $$ 2k^2-4k+1 \le 2^p $$ であれば、正当な解法と言えます。$2^p$ を固定した場合は、 $$ k \le 1+\sqrt{\frac{1+2^p}{2}} $$ であればよいですね。$k\le 2^{(p-1)/2}$ であれば十分で、Observation 7 とも合致していそうです。
Corollary 10: 任意の $(x_1, y_1, x_2, y_2)\in A_{10^9}^4$ に対して $x_1/y_1 \ne x_2/y_2$ ならば $x_1\oslash_{64} y_1 \ne x_2\oslash_{64} y_2$ が成り立つ。
Proof. $10^9\lt 2^{(64-1)/2}$ より従う。$\qed$
$2^{(64-1)/2} \approx 3.037\cdot 10^9$ なので、あまり余裕はないですね。特に、任意の符号なしの $32$-bit 整数が入力される場合はだめになっちゃいます。
Corollary 11: $i_{113} = \ceil{2^{56}+2^{27.5}} = 72057594227740468$ が成り立つ。
Proof. $i\ge 2^{(113-1)/2} = 2^{56}$ について順に探索した。$\qed$
Observation 7 でいうところの $2^{(p-3)/4}$ の項はよくわかりません。下記の記事の Observation 11 と似たような事情があるのかもしれないですし関係ないかもしれません。あんまり考えていません。
あとがき
こういう話題になると「やはり eps を使って比較するのがよい」などと安易に仰られる方がしばしばおられますが、その eps の値を決めるにあたってちゃんと議論をしましたか？というのは気になってしまいます。
また、「long double で無理やり通した」「double でやったら落ちたから long double にしたら通ったｗ」のようなツイートをしばしば見かけ、微妙な気持ちになります。
冒頭で見たように、long double でも反例がある解法もあるため、ちゃんと考えて実装することが大事ですね。
「式変形で誤差の出方が変わるなんて非本質すぎてつまらん」といった旨のツイートも見ますが、浮動小数点型においては式変形が本質であることも多く、態度を改めてほしいな〜と思います。浮動小数点型を使わされるの自体が非本質と言われると、まぁそうかもな〜という気持ちにはなります。ただ、よい精度の計算が（固定長の整数などのように）高速にできるとは思わないでほしいですね。
よく言われる「結局は Decimal が至高」みたいなのもさすがに飽きました。$10$ 進法を用いていることが本質なのではなく、Decimal のデフォルトの精度がそれなりに高いがちなことが本質で、それに気づかず $10$ 進法のおかげだと思っているのでは？と一瞬思ったのですが、別にそんなことはないかもしれません。
なんにせよ、面白い結果が得られて満足です。知っておいて役に立つかはわからないですが、役に立つかは正直どうでもいいと思っているところはありますね。 「浮動小数点型だからなんでもうまくいかないはずだ」のような思考停止をやめたい、みたいなモチベーションはあるかもしれません。 結局、大半は整数に関する議論をしているので、整数の分野が得意な人は活躍できるような気がしています。そういう人々にそういうモチベーションがある気はしないのでかなしいですね。
おわり
おわりです。
