小数点以下が高々 2 桁まである値が十進表記で与えられたとき、「それを浮動小数点数型に parse してから 100 倍する」ことで、元の値の小数点を 2 つぶん動かしたものにできますか？という話です。
>>> print(f'{float("10.00")*100:.100g}') 1000 >>> print(f'{float("10.01")*100:.100g}') 1001 >>> print(f'{float("10.1")*100:.100g}') 1010
お、いけそう？
>>> print(f'{float("10.03")*100:.100g}') 1002.9999999999998863131622783839702606201171875 >>> print(f'{float("10.05")*100:.100g}') 1005.0000000000001136868377216160297393798828125
すみません、失敗してしまいました。ということで不可能ということがわかりました。いかがでしたか？
もちろん、これで終わってはつまらないので、もう少し有益な話をします。たとえば、次のような方針たちが自然に思いつくと思います。
- 100 倍してから、ある微小値を足し、小数部分を切り捨てる
- ある微小値を足してから 100 倍し、小数部分を切り捨てる
- 100 倍した後、最も近い整数に丸める
これらの正当性について考えていきたいです。100 倍に限らず、実際に出題された例に応じた桁数を考えます。
導入
整数型以外を好まない皆さまにおかれましては、
>>> (lambda x, y: 100 * x + y)(*map(int, input().split('.'))) # 10.03 1003
や
int ipart, fpart; scanf("%d.%d", &ipart, &fpart); // 10.03 int x = ipart * 100 + fpart; // 1003
のように、「小数点で分割してから整数として parse し、ﾅﾝﾄｶ倍して足す」という方針を取ることが多いでしょう。
Python においては int('10.05'.replace('.', '')) のように小数点を消す方針でやることもありますが、桁数が固定でない場合は厄介です。
もちろんこの方法が悪いとは思いませんし、競技プログラミングにおいて初心者ができるようになるべき方針は、このような整数型を用いる方法だと思います。
藍染惣右介「整数型とは 整数型に縋らねば生きて行けぬ者の為にあるのだ」
BLEACH #47, p. 50 (407. DEICIDE9)
浮動小数点数型でやった場合でも正当性を示せる人の方がかっこいいと思うので、そうなりたいです。
前提・記法
下記を先に読むとよいかもしれません。証明などは詳しく追わなくても大丈夫だと思います。
ここで定義されている $\roundp x$, $\hfloor x$, $x\otimes y$ に関して、同様のものを用います。また、$x\oplus y \triangleq \roundp{x+y}$ も用います。
例題
- ABC 169 C: x.yy (1.2)
- ABC 226 A: xx.yyy (2.3)
- ABC 191 D: xxxxx.yyyy (5.4)
- 符号や上側の境界が他とは異なるので注意
- AGC 047 A: xxxx.yyyyyyyyy (4.9)
両方向の反例（少し小さくなるケース・少し大きくなるケース）もそれぞれ挙げておきます。
# 1.2 >>> print(f'{1.13*100:.100g}') 112.9999999999999857891452847979962825775146484375 >>> print(f'{1.09*100:.100g}') 109.0000000000000142108547152020037174224853515625 # 2.3 >>> print(f'{16.002*1000:.100g}') 16001.999999999998181010596454143524169921875 >>> print(f'{16.001*1000:.100g}') 16001.000000000001818989403545856475830078125 # 5.4 >>> print(f'{10000.0009*10000:.100g}') 100000008.99999998509883880615234375 >>> print(f'{10000.0004*10000:.100g}') 100000004.00000001490116119384765625 # 4.9 >>> print(f'{1024.000000006*1e9:.100g}') 1024000000005.9998779296875 >>> print(f'{1024.000000011*1e9:.100g}') 1024000000011.0001220703125
これらの例は binary64 の形式で roundTiesToEven による丸めを前提としています（いわゆる C++ で言う double のデフォルトの挙動）が、long double や __float128 などを使ったり、別の丸めモードにした場合でも反例は存在します*1。
なお、x.ye9 のようなリテラルの形だと、丸めが $\roundp{x.y}\otimes 10^9$ ではなく $\roundp{x.y\times 10^9}$ のように行われるため、結果が異なる場合があることに注意してください。
今回の問題設定では「文字列を一度 parse してから数倍する」であり、一旦 parse した時点で丸めが生じています。
>>> print(f'{1024.000000006*1e9:.100g}') 1024000000005.9998779296875 >>> print(f'{1024.000000006e9:.100g}') 1024000000006
考察
上で述べた方針たちについて考えていきます。
方針 1
Claim 1: $p=53$ のとき、$\delta = \roundp{0.003}$ とすると、任意の $0\le m\le 10^{13}$ に対して $$\floor{\roundp{m\times 10^{-k}} \otimes 10^k\oplus \delta} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。
Proof
$\roundp{m\times 10^{-k}} = m\times 10^{-k}+\varepsilon$ と書ける。ここで、$|\varepsilon|\le \hfloor{m\times 10^{-k}}\cdot 2^{-p}$ を満たす。
$$ \begin{aligned} \roundp{m\times 10^{-k}}\times 10^k &= (m\times 10^{-k}+\varepsilon)\times 10^k \\ &= m + 10^k\varepsilon \end{aligned} $$ であり、 $$ \begin{aligned} |10^k\varepsilon| &\le 10^k\cdot\hfloor{m\times 10^{-k}}\cdot 2^{-p} \\ &\le 10^k\cdot(m\times 10^{-k})\cdot 2^{-p} \\ &= m\cdot 2^{-p}. \end{aligned} $$ すなわち、 $$ m(1-2^{-p})\le\roundp{m\times10^{-k}}\times10^k\le m(1+2^{-p}) $$ であり、 $$ \begin{aligned} \roundp{m\times10^{-k}}\otimes10^k &\ge m(1-2^{-p})-\hfloor{m(1-2^{-p})}\cdot 2^{-p} \\ &\ge m(1-2^{-p})-m(1-2^{-p})\cdot 2^{-p} \\ &= m(1-2^{-p})^2 \\ %&= m(1-2^{1-p}+2^{-2p}) \\ &= m - m\cdot(2^{1-p}-2^{-p}). \end{aligned} $$
$m\le10^{13}$, $p=53$ とすると、 $$ \begin{aligned} \roundp{m\times10^{-k}}\otimes10^k &\ge m-10^{13}\cdot (2^{-52}-2^{-106}) \\ &\gt m-10^{13}\cdot 2^{-52} \\ &\gt m - 0.0022205 \end{aligned} $$ となる。上方向からも、同様に $$ \roundp{m\times10^{-k}}\otimes10^k \lt m + 10^{13}\cdot 2^{-52} $$ と抑えられることがわかる。
そこで、簡単のため $\delta = \roundp{0.003}$ とし、 $$\roundp{m\times10^{-k}}\otimes10^k\oplus\roundp{0.003}\lt m+1$$ を示せばよい。丸めの範囲を考えて $$\roundp{m\times10^{-k}}\otimes10^k+\roundp{0.003}\lt (m+1)\cdot (1-2^{-53})$$ を示す。 $$ \begin{aligned} &\phantom{{}={}} (m+1)\cdot(1-2^{-53}) - \left(\roundp{m\times10^{-k}}\otimes10^k\right) \\ &\gt (m+1)\cdot(1-2^{-53}) - (m + 10^{13}\cdot 2^{-52}) \\ %&= m\cdot(1-2^{-53}) + (1-2^{-53}) - m - 10^{13}\cdot 2^{-52} \\ %&= -2^{-53}m + (1-2^{-53}) - 10^{13}\cdot 2^{-52} \\ &\ge -10^{13}\cdot2^{-53} + (1-2^{-53}) - 10^{13}\cdot 2^{-52} \\ &\gt 0.9 \\ &\gt \roundp{0.003}. \quad\qed \end{aligned} $$
note: $\roundp{0.003} = \tfrac3{1000}\cdot(2^{61}+48)\times 2^{-61} \gt 0.003$ です。
方針 2
Claim 2: $p=53$ のとき、$\delta = \hfloor{\roundp{m\times 10^{-k}}}\times 2^{1-p}$ とすると、任意の $0\le m\le 10^{13}$ に対して $$\floor{(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。
なお、$\roundp{m\times 10^{-k}} + \delta$ は $\roundp{m\times 10^{-k}}$ より大きい最小の浮動小数点数であり、nextafter() や next_up() などの関数を用いることで得ることができる。
Proof
$\roundp{m\times 10^{-k}} = m\times 10^{-k}+\varepsilon$ と書ける。ここで、$|\varepsilon|\le \hfloor{m\times 10^{-k}}\cdot 2^{-p}$ を満たす。 よって、 $$ \begin{aligned} (\roundp{m\times 10^{-k}}+\delta)\times 10^k &\le (m\times 10^{-k}+\delta+\varepsilon)\times 10^k \\ &= m + (\delta+\varepsilon)\times 10^k \end{aligned} $$ となる。ここで、 $$ \begin{aligned} \delta+\varepsilon &\ge \hfloor{\roundp{m\times10^{-k}}}\times 2^{1-p} - \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\ge \hfloor{m\times10^{-k}}\times 2^{-p} - \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\ge 0 \end{aligned} $$ より、$(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k\ge m$ は成り立つ。また、 $$ \begin{aligned} \delta+\varepsilon &\le \hfloor{\roundp{m\times10^{-k}}}\times 2^{1-p} + \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\le \hfloor{m\times10^{-k}}\times 2^{2-p} + \hfloor{m\times 10^{-k}}\times 2^{-p} \\ &\le (m\times10^{-k})\times 2^{2-p} + (m\times 10^{-k})\times 2^{-p} \\ &\le 5(m\times 10^{-k})\times 2^{-p} \\ \end{aligned} $$ より、 $$ \begin{aligned} m + (\delta+\varepsilon)\times 10^k &\le m + 5m\cdot 2^{-p} \\ &\le m + 5\cdot 10^{13}\cdot 2^{-53} \\ &\lt m + 0.005552 \\ \end{aligned} $$ となる。方針 1 同様、これは $(m+1)\cdot(1-2^{-53})$ より小さいことが示せるので、$(\roundp{m\times 10^{-k}} + \delta) \otimes 10^k\lt m+1$ となる。$\qed$
方針 3
Claim 3: $p=53$ のとき、任意の $0\le m\le 10^{13}$ に対して $$\rounded{\roundp{m\times 10^{-k}}\otimes 10^k} = m$$ が成り立つ。ただし、$k$ は overflow・underflow が起きない範囲で任意の非負整数とする。
なお、実数 $x$ に対して、$x$ に最も近い整数を $\rounded x$ によって表す。round() などの関数で得られる。
ここでは、最も近い整数が複数あるケースについては考えなくてよいことが示せる。
Proof
Claim 1 の証明の中で、該当の範囲において $$ m - 0.0022205 \lt \roundp{m\times10^{-k}}\otimes10^k\lt m + 0.0022205 $$ が成り立つことを示した。これで十分である。$\qed$
方針 4
$10^k$ を掛けた後に nextafter() を行い、切り捨てる方針ではどうでしょう。証明・反証は読者への課題とします。
気が向いたら追記するかもしれません。
関連文献
あとがき
競プロの文脈では $10^{13}$ 個の全探索は絶望的ですが、こうした文脈で多少の時間を掛けてもよい状況であれば $10^{13}$ 個は全然いける寄りの個数なので、一応ちゃんと全てテストしています。特に、浮動小数点数型に関する簡単な演算であれば比較的軽い気がします。
方針 1–3 それぞれについて、1 秒あたり $10^9$ 回程度のテストができています。下記がテスト全体・ケースあたりにかかった時間です (Rust)。
- 方針 1 → 2:33:16 (0.9196 ns per case)
- 方針 2 → 4:00:15 (1.4415 ns per case)
- 方針 3 → 2:29:34 (0.8974 ns per case)
Rust に限らず C++ でも試してみましたが、.next_up() std::nextafter() は重いのかな？という感じでした。
AtCoder 環境でも $2\times 10^9$ ケース回したところ、方針 1 が 2359 ms、方針 2 が 8771 ms、方針 3 が 4760 ms (std::round()), 6893 ms (std::lround()) でした。プロセッサだったり最適化のかかり方などにもよるのかもしれません。そこの調査に関してはこの記事の範囲外という気がしています。
ともかく、こうした見積もり方がわかってくると、競プロ er あるあるであるところの
epsってよくわかんないけど、強い人が1e-6とかにしてるし、とりあえず決め打ちでテンプレに#define eps 1e-6って書いてる
のようなものが愚かに感じられるのではないでしょうか（この記事を読んで「どうやら eps は 3e-3 で決め打ちすればいいらしい」と解釈する人がいたら困ります。考えて読んでほしい）。
にぶたんの上限を常に 1e9 とか 1e18 とかで決め打ちするのが愚かだと思わないタイプの人とは思想が合わないと思うので、いちいち反論してくれなくても大丈夫です。
また、冒頭ではいくつかの成功する例も挙げましたが、浮動小数点数あるあるの「いくつかのケースではなんかうまくいく」に対する注意喚起でもあります。 貪欲法と並び、浮動小数点数型の演算も「未証明の方法は大抵誤っている」側だと思うので、証明の方針めいたものすらない状態で「いくつかテストが通ったからよさそう」というのは控えた方がよいです。
さて、一旦書こうと思っていた問題提起系の紹介記事は済ませたので、次は基礎の話を書こうかなと思います。少し期間が空くかもしれません。
これまで書いたものたち (👉 浮動小数点数):
おわり
お疲れさまでした。そろそろ Advent Calendar の時期ですね。架空 Advent Calendar 2023 の -18 日目の記事でした。
*1:入力が非負として、roundTowardPositive を使いつつ切り捨てればよいのでは？と思いつつ、それができる言語ばかりではないですね。
