はじめに
先日、接点QB氏による以下のツイートを見かけた：
??????
— 接点QB (@setten_QB) February 10, 2021
俺の脳みそがバグってるのか？？ pic.twitter.com/zyzzdb3ijS
興味深いので少し調べてみたというわけである（少々長い記事）。
numpyのarange関数の振る舞い
numpyのarange関数のステップ数に0.1を指定したときの振る舞いについて、
np.arange(1.0, 1.5, 0.1)
とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4までのarrayが得られた。array要素の最大値1.4は引数stopの1.5より小さいというのが、想定される挙動である。
一方、
np.arange(1.0, 1.6, 0.1)
とすると、0.1刻みで1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6のように、array要素の最大値が1.5を超えて1.6に達したので困惑したということである（「1.0, 1.1, 1.2, 1.3, 1.4, 1.5」かと思ったのに！）。
原因を探るため、numpyのマニュアルのarangeの項を見てみた。 numpy.org
すると、戻り値のarrayについて以下の記述が見つかった：
For floating point arguments, the length of the result is ceil((stop - start)/step).
つまり、浮動小数点数をarange関数の引数に与えた場合、array長がceil（＝小数点以下切り上げ）を伴う式で計算されるということだった。
手元で挙動を試してみると、以下のようになった。
こんな感じ？ pic.twitter.com/iJVVV91oUg
— mat (@ballforest) February 10, 2021
つまり、(1.5 - 1.0) / 0.1の計算結果が5.0となる一方で、(1.6 - 1.0) / 0.1の計算結果が6.000000000000001となり6.0にならない。
したがってceil関数を通したとき、
ceil((1.5 - 1.0) / 0.1) は5になるが、ceil((1.5 - 1.0) / 0.1) は7になった
というわけである。
問題は「なぜこのような計算結果が得られたのか」ということである。原因を一言で言えば「浮動小数点数演算に由来する計算誤差」ということに尽きるが、本記事ではこれを詳細に追ってみた（浮動小数点数の復習問題！）。
浮動小数点数の表現形式と演算
さて、浮動小数点数のIEEE754倍精度形式は次式で表される。 $$(-1)^{sgn} \times 2^{(exp - 1023)} \times (1+frac×2^{-52})$$ $sgn$は符号部、$exp$は指数部、$frac$は仮数部である（ここでは10進数表現）。
この表現形式を元にして、計算を追っていこうという方針である。 ここでは以下の記事を参考にした。 qiita.com
浮動小数点数の倍精度表現
検証にあたり、接点QB氏の具体例に基づいて、対象とする浮動小数点数を
0.1, 1.0, 1.5, 1.6
の4つに限定する。
0.1の倍精度表現
0.1の倍精度表現は
符号部 = 0 指数部 = 01111111011 仮数部 = 1001100110011001100110011001100110011001100110011010
となり、これらを10進数に変換すると
符号部 = 0 指数部 = 1019 仮数部 = 2702159776422298
となる。0.1の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1019 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{-4} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-56} \\ &= 7205759403792794 / 72057594037927936 \end{aligned} $$
1.0の倍精度表現
1.0の倍精度表現は
符号部 = 0 指数部 = 01111111111 仮数部 = 0000000000000000000000000000000000000000000000000000
であり、10進数に変換すると
符号部 = 0 指数部 = 1023 仮数部 = 0
となる。1.0の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 0 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \\ &= \frac{2^{52}}{2^{52}} \\ &= 4503599627370496 / 4503599627370496 \end{aligned} $$
1.5の倍精度表現
1.5の倍精度表現は
符号部 = 0 指数部 = 01111111111 仮数部 = 1001100110011001100110011001100110011001100110011010
となり、これらを10進数に変換すると
符号部 = 0 指数部 = 1023 仮数部 = 2251799813685248
となる。 1.5の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2251799813685248 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2251799813685248) \times 2^{-52} \\ &= (4503599627370496 + 2251799813685248) \times 2^{-52} \\ &= 6755399441055744 / 4503599627370496 \end{aligned} $$
1.6の倍精度表現
1.6の倍精度表現は
符号部 = 0 指数部 = 01111111111 仮数部 = 1001100110011001100110011001100110011001100110011010
となり、これらを10進数に変換すると
符号部 = 0 指数部 = 1023 仮数部 = 2702159776422298
となる。 1.6の倍精度表現は以下のように計算できる。 $$ (-1)^{0} \times 2^{1023 - 1023} \times (1 + 2702159776422298 \times 2^{-52}) $$ $$ \begin{aligned} &= 1 \times 2^{0} \times (2^{52} + 2702159776422298) \times 2^{-52} \\ &= (4503599627370496 + 2702159776422298) \times 2^{-52} \\ &= 7205759403792794 / 4503599627370496 \end{aligned} $$
倍精度表現に基づく浮動小数点数演算
以上で倍精度表現が揃ったので、検証に入ることができる（前置き長すぎワロタ）。 表に示しておく。
| 浮動小数点数 | 倍精度表現 |
|---|---|
| 0.1 | 7205759403792794 / 72057594037927936 |
| 1.0 | 1 (=4503599627370496 / 4503599627370496) |
| 1.5 | 6755399441055744 / 4503599627370496 |
| 1.6 | 7205759403792794 / 4503599627370496 |
以上の表に基づいて、(1.5 - 1.0) / 0.1 を計算すると
((6755399441055744 - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936) = 5
が得られた。ぴったし5である。
一方、(1.6 - 1.0) / 0.1 については
((7205759403792794 - 4503599627370496) / 4503599627370496) / (7205759403792794 / 72057594037927936) = 6.000000000000001
と計算できた。
以上の計算結果に基づいてceilは小数点以下を切り上げるため、冒頭で示した挙動が説明できたわけである。
まとめ
本記事ではarange関数の引数のひとつであるstep数に浮動小数点数を適用した場合に、 戻り値のarrayの要素数がどのように変わりうるかを、具体的な数値例を元に調べてみた。
Numpyのマニュアルによれば、arangeによる戻り値のarrayの要素数はceil((stop - start)/step)で計算され、これらstart, stop, stepはarange関数の引数として与えられる。
接点QB氏のツイートをもとに、arangeの引数が具体的に
- (1) start=1.0, stop=1.5, step=0.1
- (2) start=1.0, stop=1.6, step=0.1
の場合に、浮動小数点数の倍精度表現に基づいて要素数を計算し、arange関数の振る舞いを検証した。 検証の結果、確かにPythonの計算結果は正しかったことが確かめられた。
Numpyマニュアルにも浮動小数点数が引数の場合はlinspace使ったほうがいいよ的な記述もあり、arangeを使うときには思い込みによるバグを埋め込まないように気をつけたい。
