pythonの除算結果が浮動小数点数になったので、必要な精度が失われた話

(注:pythonの話を書いていますが、浮動小数点の精度の話は他のプログラミング言語でも同様だと思います。 pythonの記法では、浮動小数点の除算演算子/であり、整数除算の演算子//であることに注意してください)

特に、大きな数の計算で答えの整数を正確に求める必要がある場合に、この問題にハマりやすい。

Pythonでは除算の結果がfloatになる

公式リファレンスより。

整数の除算結果は浮動小数点になりますが、整数の切り捨て除算結果は整数になります
https://docs.python.org/ja/3/reference/expressions.html#binary-arithmetic-operations

公式チュートリアルより。

除算 (/) は常に浮動小数点数を返します。 https://docs.python.org/ja/3/tutorial/introduction.html#numbers

整数の除算結果は浮動小数点、float型になる。実際に割り切れるかどうかとは無関係である。

>>> 9/3
3.0
>>> 10/3
3.3333333333333335
>>> type(9/3)
<class 'float'>
>>> type(10/3)
<class 'float'>

floatが持っている精度が、必要としている精度よりも小さいと、数値が不正確になる

浮動小数点数はその性質上、数を正確に表現できない。 例えば、上で計算した10/3を正確に表現しようとしたら無限の桁が必要になってしまうので、有限桁で数値を丸める必要がある。 当然ながら、一般的にはこれによって問題は発生しない。

しかし、多くの桁数に対して精度が必要な場合は別である。

浮動小数点型はたいていは C の double を使って実装されています https://docs.python.org/ja/3/library/stdtypes.html#numeric-types-int-float-complex

というわけでIEEE 754 - Wikipediaを見てみよう。(実際のところC言語のdoubleが必ずIEEE 754に準拠するかは、まだ調べていない……)
倍精度浮動小数点では、64ビットを使って数を表現し、その内訳は次のとおりである。

  • 符号ビット (sign): 1 ビット
  • 指数部 (exponent): 11 ビット
  • 仮数部 (fraction): 52 ビット

精度を求める上で重要なのは仮数部のビット数 = 52ビットである。これに加えて、「最上位ビットは1に決まってるんだから書かなくてもいいでしょ、省略するわ」という省略された1ビット(ケチ表現)がある。合計53ビットがdouble型の精度となり、10進に直した精度は53 * log10 2 =15.95桁となる。

つまり、float型の精度は約16桁である。これより高い精度は表現できない。

桁数が16桁よりも多い割り算で、floatとintの場合に何が起きるか、例を見てみよう。

>>> 111111111111111111 // 9
12345679012345679
>>> 111111111111111111 / 9
1.234567901234568e+16
>>> int(111111111111111111 / 9)
12345679012345680

1を18個並べた数を9で割ると、ちょうど割り切れて商は12345679012345679である。これは17桁の数である。
切り捨て除算の演算子を使うと、計算結果はintのままなので正確な値となる。
しかし、(不必要に)浮動小数点の除算をすると、計算結果はfloatになるのでそこでの精度は10進数で約16桁である。17桁よりも少ないので、計算結果を必要な精度で表現できない。 その後でint型にキャストしたとしても、失われた精度は戻らない。直接割った値とは異なる値になってしまう。
17桁の数を正確に求めようと思ったら、浮動小数点の除算を経由してはいけないのだ。

整数だけで済むなら、整数で閉じるように演算を工夫する

一般に、最終的な結果が整数であれば、上手く計算して途中結果がfloatにならないように、intだけになるようにすればよい。

実のところ、先日の自分がつまづいたのは、最終的な結果が整数のケースだった。 整数a,bに対して、aをbで割った数を切り上げた数を求めようとしていたのだ。

これを実行するには2つの方法がある。
浮動小数点数を経由してよければ、 math.ceil(a/b) とすればよい。
また整数除算だけで実行するには、(a+b-1)//bとすればよい。

計算速度の観点を除外して結果の正しさだけを考えた場合、小さな数ならば両者に違いは無い。
したがって、どちらを使っても正しい結果が得られる。

>>> import math
>>> a,b = 10, 2
>>> math.ceil(a/b)
5
>>> (a+b-1)//b
5
>>> math.ceil(a/b) == (a+b-1)//b
True
>>>

しかし、桁数が増えるとmath.ceil(a/b)の値は不正確になる。除算によっていったんfloatになったときに、必要な精度が失われてしまうからだ。
したがって、正しい結果を得るためには(a+b-1)//bのほうを使わねばならない。
下の例では除算の結果が10進17桁なので、a/bの結果が不正確になっている。

>>> a,b = 98765432101234567, 2
>>> math.ceil(a/b)
49382716050617280
>>> (a+b-1)//b
49382716050617284
>>> math.ceil(a/b) == (a+b-1)//b
False

より高い精度が欲しいならdecimalモジュール

らしいですが、こっちはまだ全く調べていません。


今日のまとめは以下の通り。

  • Pythonでは除算の結果がfloatになる
  • floatの精度が、必要としている精度よりも小さいと、期待する結果が得られない
  • 途中の結果も整数だけで済むなら、整数で閉じるように演算を工夫する

2020年6月28日 追記

競技プログラミングの問題を解いていてこの問題にぶち当たった、という話を書いていなかった。 何で競技プログラミングに言及しなかったんだろ。

最初にこの問題に引っかかったのは 、ARC062の「AtCoDeerくんと選挙速報」である。
C - AtCoDeer and Election Report
「は!? 何でWAになるんだ、合ってるだろ」と間違いの原因が分からず、テストケースを見に行ってやっと原因が分かった。

結構な人がこの問題では上記のバグを踏んでハマっている。
俺が一言もこの問題に言及しなかったのに、この記事にリンクを張ってくれてる解説記事がある。
ARC-062 C-AtCoDeerくんと選挙速報: ceil関数の割り算に気をつけよう - Qiita

また、Pythonに限らず、C++でも同じ罠があることが分かった。(PythonでもC++でも浮動小数点の仕様はほぼ同じだから、同様の結果になる)。
ARC 062 C - AtCoDeerくんと選挙速報 / AtCoDeer and Election Report - ykmakuのブログ
のceilを使ってWAになっている部分である。

また、AGC044の「Pay to Win」でも同様の問題になる。
A - Pay to Win
解説ではfloorとceilを使っているが、その通りに実装してはいけない。
Nが1018以下、という制約である。floorとceil関数を使うと、この問題に引っかかる。

わざとfloor, ceil関数を使って書いてみると、下記のようにWAになります。
Submission #14813496 - AtCoder Grand Contest 044

正解はこちら。コード実行時間も結構変わってくるのね。
Submission #14823409 - AtCoder Grand Contest 044