いものやま。

雑多な知識の寄せ集め

Pythonの四捨五入で四苦八苦した話。(その3)

前回のあらすじ:Rubyの四捨五入すごい。PyPIにあったパッケージどれもダメ。

Rubyでの四捨五入のアルゴリズム

さて、そんなわけで、Pythonがどうしようもなかったので、なぜか誤差なく四捨五入できてたRubyでの実装を調べてみた。

Rubyでの四捨五入は、numeric.cのflo_round()で実装されている。

flo_round()ではいくつかの難しいケースへの対処も入っていて、通常処理の本体はマクロROUND_CALL()から呼び出されるround_half_up()。 これを解読してPythonで書き直すと、次のようになる:

def myround3(number, ndigits=0):
    shift_amount = math.pow(10, ndigits)

    if number >= 0:
        shifted = math.floor(number * shift_amount)
        bound = (shifted + 0.5) / shift_amount
        if number >= bound:
            shifted += 1
    else:
        shifted = math.ceil(number * shift_amount)
        bound = (shifted - 0.5) / shift_amount
        if number <= bound:
            shifted -= 1

    return shifted / shift_amount

これを実際に試してみると、次のような感じ:

myround3(0.5)  # => 1.0, OK
myround3(1.5)  # => 2.0, OK
myround3(-0.5)  # => -1.0, OK
myround3(-1.5)  # => -2.0, OK
myround3(12345.67, -1)  # => 12350.0, OK
myround3(-12345.67, -1)  # => -12350.0, OK
myround3(1.255, 2)  # => 1.26, OK
myround3(-1.255, 2)  # => -1.26, OK

基本的な四捨五入はもちろん、鬼門だった1.255を小数第三位で四捨五入するのも正しくできてるのが分かる。

では、なぜこのアルゴリズムだとうまくいくのか?

ポイントは、シフト前の状態で判断しているということ。

これまでのアルゴリズムだと、1.25xxxをシフトして125.xxxにし、そこに0.5を足して126を超えるかどうかで切り上げ/切り捨てを判断してた。 つまり、シフト後の状態で判断をしていたことになる。 けど、これだとシフトで誤差が露呈してしまい(1.255をシフトすると125.499...になる)、うまくいかなかった。

一方で、Rubyでのアルゴリズムは1.25xxxを四捨五入するときの境界値(この場合は1.255)を作り出し、これとの比較をすることで切り上げ/切り捨てを判断している。 つまり、シフト前の状態のままで判断をおこなっている。 こうすると何が嬉しいかというと、比較のときに元の値をそのまま使えるので、計算での誤差がそこに乗らない

もちろん、境界値を作る計算では誤差が乗るんだけど、面白いのは、誤差が乗ったとしてもその誤差は元の数値と同程度にしかならないということ。 この非対称性が面白い。

これは具体的に見てみるとよく分かる:

1.255 * 100  # => 125.49999999999999, 誤差が露呈する
125.5 / 100  # => 1.255, 誤差が乗ってない(ように見える;内部では1.255と同じ誤差が乗ってる)

このように、シフトには誤差が乗るけど、アンシフトには誤差が乗らないという性質をうまく使ったのがRubyでの実装となっていた。 いやー、賢い・・・

さらに厳しいケースへ

このRubyリポジトリでは単体テストも用意されてて、RSpecで書かれたspec/ruby/core/float/round_spec.rbおよびtest-unitで書かれたtest/ruby/test_float.rb, test/ruby/test_numeric.rbなどがある。

これらの中にはバグチケットの番号とともに書かれた厳しいテストケースなどもある。

たとえば、「42.0を小数点以下308桁になるように四捨五入したら42.0になること」みたいな。 これをさっきのmyround3()でやってみると、こうなる:

myround3(42.0, 308)  # => OverflowError: cannot convert float infinity to integer

これは、10の308乗もの数を42.0に掛けるので、浮動小数点数で表現できる桁を超えてしまって(オーバーフロー)、math.floor()が切り捨てできないよってなるのが原因。 こうなる場合、四捨五入する桁はずっと小さい桁になってるから元の数値に影響を与えないことになるので、元の数値をそのまま返すことになる。 そして、オーバーフローが発生するかどうかは、ざっくり次のようなロジックで書ける:

def will_overflow(number, ndigits):
    _, num_binexp = math.frexp(number)
    num_exp = num_binexp / 3.322259  # log_2(10) = 3.322259
    return (num_exp + ndigits) >= sys.float_info.dig

sys.float_info.dig浮動小数点数を10進数で表現したときに信頼できる桁数を表す定数。 つまり、この桁数を Kとすると、10^K以上の数は表現できないことになる。 そして、元の数値が 10^mと表現でき、シフトする量が 10^nなのだとすれば、これを掛けた値は 10^m \times 10^n = 10^{m+n}となる。 そこで、 10^{m+n} \ge 10^K、すなわち m + n \ge Kとなったら、オーバーフローという感じ。 (なお、math.frexpでは2進数での桁数が取れるので、それを10進数の桁数に直す計算が最初に入ってる)

逆に、「42.0を230桁で四捨五入したら0.0になること」というのも。 これをさっきのmyround3()でやってみると、こうなる:

myround3(42.0, -(2**30))  # => ZeroDivisionError: float division by zero

これは逆に浮動小数点数で表現できる一番小さい幅よりも小さい値になってしまい(アンダーフロー)、値が0に落ちてしまってゼロ除算を発生させてしまったのが原因。 これも同様にアンダーフローが発生するか((num_exp + ndigits) < 0か)を確認し、発生していたら、それは極端に大きな桁で四捨五入していることになるので、0を返せばいい。

さらに、さきほどのsys.float_info.digに関連して、これより大きな桁については信頼できる数字ではなくなってしまうため、次のようなことが起こる:

for i in range(1, 32):
    print(f"[{i:02}] {myround3(0.5, i)}")
# 以下の出力になる
# [01] 0.5
# ...(省略)
# [15] 0.5
# [16] 0.5000000000000001
# [17] 0.5
# ...(省略)
# [31] 0.5

実行したマシンではsys.float_info.digが15なんだけど、それを超えた16桁のときに誤差が入り込んでしまってる。

こうなるともう浮動小数点数では対処しきれなくて、Rubyでは有理数を使った計算に切り替えて対処している。 Pythonの場合は10進数を10進数のまま扱えるようにするdecimalを使うといい。

パッケージへ

とまぁそんな感じで、たかが四捨五入、されど四捨五入。 ちゃんとやろうとすると、かなり大変なことが分かる。

なので、パッケージを作って単体テストRubyから移植し、マトモで堅牢な四捨五入を使えるようにした。 それがsaneroundパッケージ。

四捨五入で四苦八苦しないためにも、ぜひ使ってほしい。

今日はここまで!

Pythonの四捨五入で四苦八苦した話。(その2)

前回のあらすじ: Pythonのビルトインのround()は普通の四捨五入とは違う偶数丸めというものだったので、floor()を使って普通の四捨五入を実装したよ。

誤差との戦い

と、まぁ普通はここまでなんだけど、いろいろ試すと「あれ?」ってなるケースが出てくる。

具体的には、1.255を四捨五入して小数第二位までにしようとしたとき。 普通に考えれば小数第三位の数字が5なので切り上げられて1.26になると思うけど、残念ながらそうはならない。

myround2(1.255, 2)  # => 1.25

ちなみに偶数丸めするround()ならこれは1.26になるはずなんだけど、試したらこうだった:

round(1.255, 2)  # => 1.25

Pythonくんさぁ・・・

それはともかく、どうしてこうなるのかというと、誤差の関係。 浮動小数点数だと10進数の数字は正確にもてない場合があって、その場合に計算するといろいろ誤差が出てしまう。

これは実際に1.255 * 100を計算してみるとよく分かる:

1.255 * 100  # => 125.49999999999999

見てのとおり、きっちり125.5になっていなくて125.499999...となってしまっている。 このせいで0.5を足しても126に上がってくれなくて、結果1.26ではなく1.25となっていた。

まぁ、誤差が原因だとしょうがない感じはある。 きっと他のプログラミング言語でも同様の制約は出るはずだし。

と思って試しにRubyでやってみた結果がこれ:

1.255.round(2)  # => 1.26

えっ、ちゃんと四捨五入できてる? なんで???

ちなみに、もちろんRubyでも内部では誤差を持ってて、1.255 * 100を計算してみるとこうなる:

1.255 * 100  # => 125.49999999999999

じゃあ、Rubyは一体どうやって四捨五入を行なっているのか。 それについてはまたあとで言及したい。

PyPIのパッケージ

まぁ何はともあれ、四捨五入が思ったよりも大変そうということが分かったので、だったら誰かがパッケージを作って公開してるだろうとPyPIを検索。 すると、案の定いくつかのパッケージが見つかった:

最初に言っておくと、これらのパッケージはダメダメなので使ったらダメ。 (宣伝になってしまうけど、これらを使うくらいなら拙作のsaneroundを使って・・・)

実際、試してみた結果がこれ:

まずはround2。

from round2 import round2

round2(0.5)  # => 1, OK
round2(1.5)  # => 2, OK
round2(-0.5)  # => -1, OK
round2(-1.5)  # => -2, OK
round2(15, -1)  # => nan, NG
round2(25, -1)  # => nan, NG
round2(1.255, 2)  # => 1.25, NG

見てのとおり、四捨五入の桁で負の数は未サポート。 また、1.255の小数第三位を四捨五入しても1.26にはならずに1.25になってしまってる。

次にmath-round。

from math_round import mround

mround(0.5)  # => 1, OK
mround(1.5)  # => 2, OK
mround(-0.5)  # => 0, NG
mround(-1.5)  # => -1, NG
mround(-2.0)  # => -1, NG(根本的な問題がありそう)
mround(15, -1)  # => 20.0, OK
mround(25, -1)  # => 30.0, OK
mround(1.255, 2)  # => 1.25, NG

こちらは負の数を四捨五入したときの挙動がおかしい。 単に絶対値で扱ってないことによる問題かと思いきや、-2を四捨五入すると-1になってるあたり、根本的な問題がありそう。 そして当然1.255の四捨五入結果は1.25。

最後にmath-round-af。

from math_round_af import get_rounded_number

get_rounded_number(0.5)  # => 1.0, OK
get_rounded_number(1.5)  # => 2.0, OK
get_rounded_number(-0.5)  # => -1.0, OK
get_rounded_number(-1.5)  # => -2.0, OK
get_rounded_number(15, -1)  # => ValueError, NG
get_rounded_number(25, -1)  # => ValueError, NG
get_rounded_number(1.255, 2)  # => 1.25, NG

こちらも四捨五入の桁で負の数は未サポートとなっている。 また、1.255の四捨五入結果はやっぱり1.25。

とまぁこんな感じで、精度の問題どころか、基本的な実装すらできてないお粗末なパッケージしかなくて、本当にダメダメ。 3つも引っ張ってくりゃ普通は1つくらいまともなパッケージありそうなもんだけど。。。

もちろん、ちゃんと探せばもしかしたらまともなパッケージがまだあるかもしれない。 でも、これはもう自分でなんとかした方が早そうだなとなった。

ということで、次はRubyの実装を調査していくことになる。

今日はここまで!

Pythonの四捨五入で四苦八苦した話。(その1)

Pythonの四捨五入は、実は普通の四捨五入になってない。

なので、普通の四捨五入をするにはちょっとコードを書く必要があったりする。 ただ、その過程でかなり苦労があったので、同じような苦労をする人がいないようにパッケージにまとめてPyPIに上げておいた。

使い方は簡単で、pip install saneroundでインストールしたあと、saneroundをインポートし(srとしておくと楽)、round()を呼び出すだけ:

import saneround as sr
sr.round(0.5)  # => 1.0

もちろん、負の数を指定してもOK:

sr.round(-0.5)  # => -1.0

第2引数に桁を指定すると、その桁になるように四捨五入する:

sr.round(1.255, 2)  # => 1.26
sr.round(1234.567, -2)  # => 1200.0

さて、結論としてはもうこのsaneroundを使ってねなんだけど、Pythonで四捨五入をしようとするとどういった苦労があるのかを以下では書いておきたい。

ビルトインのround()関数

Pythonだって一応はプログラミング言語の端くれなので、四捨五入をするround()関数くらいは標準で入ってる。 ということで、四捨五入しようと思ってこのround()関数に手を出してしまうと、大火傷する。

何が起こるのかは0.5を四捨五入してみれば一目瞭然で、0.5を四捨五入すれば当然1になるとみんな思うでしょ?

はい、

round(0.5)  # => 0

0ですね。

「あー、そういうことね、0.5未満じゃなくて0.5以下が下に丸められるのね、理解」

round(1.5)  # => 2

はぁ? 頭おかしいの???

・・・まぁ、Pythonが頭おかしいのはもう仕方ないとして、この仕様を知らずに使ってたりすると、大問題になる危険性がある。 そこだけは注意する必要がある。

なお、この仕様自体は「偶数丸め」と呼ばれるもので、境界にあたる数(0.5や1.5など)を四捨五入したときには、丸られた結果の最小桁が偶数になるように切り上げ/切り捨てするというもの。

たとえば、以下のような感じ:

round(0.5)  # => 0, 一の位が偶数
round(1.5)  # => 2, 一の位が偶数
round(1.25, 1)  # => 1.2, 小数第一位が偶数
round(1.35, 1)  # => 1.4, 小数第一位が偶数
round(25, -1)  # => 20, 十の位が偶数
round(35, -1)  # => 40, 十の位が偶数

なんでこんな変な仕様を標準にするかなぁと思うんだけど、まぁ仕方ない。 だってPythonだし。。。

「普通の四捨五入」の自作

ということで、普通の四捨五入をしようと思ったら、普通は自分でコードを書くことになる。 四捨五入はfloor()関数使えば実装できることはよく知られてるし。

たとえば、よくある単純な実装は以下のようなもの:

import math

def myround(number, ndigits=0):
    shift_amount = 10 ** ndigits
    shifted = number * shift_amount
    return math.floor(shifted + 0.5) / shift_amount

やっているのことは、

  1. 四捨五入したい桁が小数第一位にくるようにする
  2. 0.5を足すことで、小数第一位が0.5以上なら一の位が1増え、それ以外なら増えないようにする
  3. math.floor()で小数部分を切り捨て
  4. 元の桁に戻して返す

これで実際に動かしてみると、以下のようになる:

myround(0.4)  # => 0.0
myround(0.5)  # => 1.0, round()だと0.0
myround(1.4)  # => 1.0
myround(1.5)  # => 2.0
myround(1.14, 1)  # => 1.1
myround(1.15, 1)  # => 1.2, round()だとなぜか1.1(バグっぽい)
myround(1.24, 1)  # => 1.2
myround(1.25, 1)  # => 1.3, round()だと1.2
myround(14, -1)  # => 10.0
myround(15, -1)  # => 20.0
myround(24, -1)  # => 20.0
myround(25, -1)  # => 30.0, round()だと20

ちゃんと四捨五入できてて、偶数丸めもされてないじゃん、めでたしめでたし。 とは残念ながらならなくて、負の数を指定するとおかしくなる:

myround(-0.4)  # => 0.0
myround(-0.5)  # => 0.0, 理想は-1.0
myround(-0.6)  # => -1.0
myround(-1.4)  # => -1.0
myround(-1.5)  # => -1.0, 理想は-2.0
myround(-1.6)  # => -2.0

これは絶対値で見ていないのが原因。 なので、次のように直してやるといい:

def myround2(number, ndigits=0):
    sign = 1 if number >= 0 else -1
    shift_amount = 10 ** ndigits
    abs_shifted = sign * number * shift_amount
    return sign * math.floor(abs_shifted + 0.5) / shift_amount

符号を確認して、絶対値に直した上で四捨五入し、元の符号に戻してやる感じ。

これで動かすとこうなる:

myround2(0.4)  # => 0.0
myround2(0.5)  # => 1.0
myround2(-0.4)  # => 0.0
myround2(-0.5)  # => -1.0

完璧。


・・・で終わってくれればよかったんだけど、そうはならないのが残念なところ。 長くなってきたので続きは次回で。

今日はここまで!

『withのススメ』というLTをしてみた。

5/25(水)に、株式会社ラクスさんの主催する勉強会「Python Tips LT会 - vol.3」でLT登壇してみた。

話したのは『withのススメ』というLT。

まぁ、Pythonの話をしつつ、ホントにしたいのはRubyはいいぞっていう話だったりするけどw

ちなみにwith文はもっと使われた方がいいと思っていて、たとえば他の発表でPythonでシリアル通信やった話があったけど、ここでもclose()の呼び出しは必須だし、あるいはセマフォで資源数の管理をした話についてもwith文で管理した方が安全といえる(資料だと正常処理とexceptの両方にrelease()があって、まだfinallyの方がマシだけど、withならacquire()/release()を書く必要すらなくなる)。 ユーザが自前でwith文を使えるようにするのはちょっと大変だから、ライブラリ側でちゃんと用意してくれるといいんだけどね。

今日はここまで!

オブジェクト指向と関数型の話。

「たまにはLTやろうぜ」というイベントがあったので参加してLTしてきた。

スライドはこちら:

なお、前回の記事(『オブジェクト脳への一歩』というLTをしてみた。 - いものやま。)で、「そもそも関数型と対比されるべきは手続き型で、オブジェクト指向ではない」ということを書いたけど、副作用への対処という部分に注目すると、オブジェクト指向と関数型ではまったく逆のアプローチをとっていると思っていて、その部分については比較して考えることもできる。 まぁ、結論としては「関数型はオブジェクト指向のサブシステムにしかなりえない」ということで、やはりオブジェクト指向が上位にある考え方だとなるんだけど。

ちなみに、なんでこの話をしようと思ったのかというと、『ミノ駆動本』こと『良いコード/悪いコードで学ぶ設計入門』を読んでて、関数型に毒されすぎだよなと思ったから。

この本、今めっちゃ売れてるらしく、自分も読んだんだけど、とにかくイミュータブル推しで、でもそれってオブジェクト指向的な考え方とはちょっとズレてる。 もちろん、変わらないものは不変であるべきなんだけど、変わるものまで不変にしてしまうのは、ちょっとやりすぎ。。。 それはデータの変更に対する責務を外部に押し付けてるだけで、外側はツラいんだよなぁ。 関数型はそうするしかないからそうしてるのであって、オブジェクト指向ではオブジェクト指向らしい対処法があるのだから、ちょっと考えてほしいところ。 この本がやたら売れてるだけに、なんでもかんでもイミュータブルの民が生まれなければいいんだけど。。。

今日はここまで!

『オブジェクト脳への一歩』というLTをしてみた。

株式会社ラクスさんの主催する勉強会で「リーダブルコード LT会 - vol.3」というのがあったので、LT登壇してみた。

話したのは『オブジェクト脳への一歩』というLT。

オブジェクト指向の本はいろいろあるけど、一番の入り口である「処理から考え始めるんじゃなくて、データから考え始める」という部分に自覚的な本って見た記憶がないよなというのがモチベーション。 といっても、自分もこれに気がついたのは最近で、意識することなくいつのまにかオブジェクト指向的な考え方に慣れていたというのが実際のところ。 でも、ここに自覚的になれると、オブジェクト指向でプログラミングするってこういうことだったんだと一気に視界が開けるんじゃないかなと思う。

オブジェクト指向の本の話

ちなみに、「オブジェクト脳」という言葉を使ってる本といえば、コレ:

評判悪くない本なんだけど、自分が読んだときの感想は、たぶんこれ読んでもオブジェクト指向できるようにはならないよなぁというもの。

あと、オブジェクト指向に関する有名な本といえば、コレ:

この本も割と残念というか、別にオブジェクト指向に限った話じゃないことをさもオブジェクト指向のおかげであるかのように話してたり、視野の狭さが気になった。

いずれにしても、処理を中心に考えるんじゃなくて、データを中心に考えるんだという発想の転換が示されてないので、オブジェクト指向的な考え方に至らないんよね。

良さげな本で頭に浮かんだのはオライリーのHead Firstシリーズの『オブジェクト指向分析設計』だけど、今見返すと、うーん、ちょっと余分な情報が多くて読みづらいかも。絶版だし、Javaのコードも古くなっちゃってるからなぁ。

ただ、ユースケースをどう書いたらいいのかとか、そこからどう設計に落とし込んでいくかなど、学べることは多いと思う。 同じHead Firstシリーズの『デザインパターン』の本も悪くない(ただしこれも絶版だし古い)。

なんかいい本があるといいんだけど。 (自分で書くしかない?)

関数型とオブジェクト指向

ちなみに、最近はオブジェクト指向なんか古くて今は関数型だよ、という人がいるけど、それは全然オブジェクト指向を理解できてない。

関数型と対比されるべきなのは手続き型で、これらは計算モデルが違っていて、関数型はラムダ計算を基礎に持ち、一方で手続き型はチューリングマシンによる計算を基礎に持つ。

対するオブジェクト指向というのは、関数型、手続き型の両方の上位にある考え方で、関数型とぶつかるものじゃない。

実際、関数型言語であるErlangの開発者の一人であるアームストロング氏は、最初オブジェクト指向を批判したものの、のちに「Erlangこそオブジェクト指向である」と言っていたりする。

結局、いい設計になってるかどうかが肝心で、オブジェクト指向言語で書いていてもオブジェクト指向的でないゴミのような設計になってれば、当然分かりにくい。 それをオブジェクト指向のせいにしないでほしいなぁ。

今日はここまで!

簡単な組版ツールを作ってみた。

技術同人誌を書こうとすると必要となるのが組版。 これまでは \TeXを使ってきたけど、いろいろとツラいので、もっと使いやすい組版システムが欲しいなと。

そんな想いを書いたのが技書博5で出した『TeXグッバイしたい本 vol.0』というコピー本。

この本自身、 \TeXは使わずにRubyで出力はしてるんだけど、酷い組版で。。。

そこからもうちょっと進歩して、今回の技術書典12で出したのが『TeXグッバイしたい本 vol.1 フォントのはなし』。

この本を書くために、前回のツールをもう少し発展させて、簡単な組版ツールを作ってみた。 まぁ、これもまだまだちゃんとした組版と呼べるものにはなってないんだけど。

そして、せっかくなのでこのツール部分を公開することにした。

上記の本ではこのコードのうちフォントを扱ってる部分について解説してる。 興味のある人は書籍の方もぜひ。

ちなみに、PDF出力してる部分や組版している部分についても続刊で説明していく予定。 いろいろ機能追加もしていきたいので、乞うご期待。

なお、リファクタリングとか設計変更でコードもガラッと変わっていく可能性が高いので、使うのは自由だけど、急に使えなくなったとか出てきてもそこは自己責任で。 将来的にはこのツールをそのまま使うのではなく、ここでのフィードバックをベースにgemを作っていきたいと思ってる。

今日はここまで!