Pythonと1/9801

先日、こんなツイートを見ました

1から100の自然数を表示させる方法というものですが、その中の一つに1/9801を小数で表示すると0.0001020304050607... というふうに2桁の整数が順番に表れるというものがありました。え、すげー！

その後しばらくしてふと、そういえばこれ、数年前にコーディング面接って文脈で見たことあったなと思い出しました。*1 コーディング面接中にこんな回答したら面接官の人も驚くだろうなあ～～～～～

というわけで、早速

最近家に籠ってばかりで気分もなかなか上がらないので、天才プログラマごっこをして遊ぶことにしました。

面 「それでは、1から100までの自然数を出力するコードを書いていただけますか？言語は自由で結構です。」
僕 「はい、ではpython3でやらせていただきます。」*2

print(1/9801)

>> 0.00010203040506070809

おっと？

面「もう少し出せますか？」
僕「はい！」

焦ることはありません。formatを利用して表示する桁数を増やしてあげましょう。なにしろ自分はこの後にも表示が続くことを知っていますから。手始めに50桁くらい出しますかね。

print("{:.50f}".format(1/9801))

>> 0.00010203040506070808982747033510918299725744873285

おっとっと？

面「面接は以上です。」

はい

というわけで、なんでも計算してくれる計算機くんも厳密な値を返してくれるわけではないことがわかりました。
通常、計算機の実装の都合上、数値は2のべき乗の和をずらしたものでしか表現できません。(小数なら 1/2, 1/4, 1/8, 1/16 ... の和。詳しくは「浮動小数」で。)

floatの仮数部は53ビットらしいので、253 ≒ 9 * 1015 より、先頭の0を除いて16桁目くらいから崩れる可能性があるということですね。 pythonでは公式ドキュメントでこの現象に付随するtipsがいくつか紹介されています。

15. 浮動小数点演算、その問題と制限 — Python 3.8.2 ドキュメント

残念。

とはいえ

脳内面接シミュレーション第1回はお祈りを食らってしまいましたが、せっかくならちゃんと最後までやりたいと思うわけです。

実は上記のドキュメントに解決方法が書いてあり。Decimalというモジュールを使うことができます。
Decimalは10進表記に寄り添った形で演算を行っており、小数点以下の計算も正確に行うことができます。

解決策を得たところで面接シミュレーション第2回を始めます。めげない心が大事です。 もう精度に気を使う必要はありませんから、いきなり50桁出しても大丈夫でしょう。

from decimal import *
a = Decimal(1)/Decimal(9801)
print("decimal: {:.100f}".format(a))

> 0.00010203040506070809101112131420000000000000000000

おやおや

実はDecimalのデフォルトでは28桁でやめることになっています。割り算の筆算を途中で四捨五入してやめたのと同じ感じでしょうか。
要求する精度はこちらで指定することができますので、あわてる必要はありません。

getcontext().prec = 50

b = Decimal(1)/Decimal(9801)
print("decimal: {:.50f}".format(b))

> 0.00010203040506070809101112131415161718192021222324

すばらしい！

いやー美しい計算結果ですね。精度は100桁でも500桁でも1000桁でも指定できます。当初の予定通り100までの自然数が欲しければ200桁ちょいを指定してあげれば解決ですね、めでたしめでたし！*3

いかがでしたか？いやあ、冒頭の状態で現実の面接を受けていたら大惨事になるところでしたね。やっぱりシミュレーションを実際におくってのはとっても大事なんですね。みなさんにもわかっていただけましたか？

おまけ1: fractions

ここからは、正確な計算結果を得るもう一つの方法であるfractionsについてです。
fractionsは小数を有理数の状態で扱うことで正確な計算を行うというモジュールです。

from fractions import *
print(Fraction(1, 13) + Fraction(1, 13)) # 1/13 + 1/13

> 2/13

ただし、小数を入力とするときは結局一度Decimalをかませるなどして「正確な0.1」を準備する必要があります。

from decimal import * 
from fractions import *

print("{:.30f}".format(0.1)) # <- 実は正確な"0.1"じゃない
print(Fraction(0.1))　　　# <- なので、これも想定した結果と異なる

print(Fraction(Decimal("0.1"))) # <-正確な"0.1" が入力になっている

> 0.100000000000000005551115123126
> 3602879701896397/36028797018963968
> 1/10

余談ですが、fractionsにはひとつだけなんとなく浮いてるメソッドが存在します。AtCoderのpython勢ならご存知かと思いますが、そうです、fractions.gcd() です。*4

前まではfractionsってmathの前身かなんかなのかな、と思っていたのですが、既約分数を得るためのメソッドを借りていただけだったのですね。

おまけ2: 浮動小数でも頑張りたい

ここまでで正確に計算する方法はわかったのですが、ちょっとだけ面倒な回り道してるという感覚もあります。

公式ドキュメントでは「もっと詳しくやりたかったらnumpyも参照してみてね」とも言っていました。

そこで、200桁は無理でももう少し精度の高い型が提供されているのでは？と思い調べてみました。が、

import numpy as np

a = np.half(1)/np.half(9801)
print("numpy half      : {:.50f}".format(a))

b = np.float(1)/np.float(9801)
print("numpy float     : {:.50f}".format(b))

c = np.double(1)/np.double(9801)
print("numpy double    : {:.50f}".format(c))

print("========================================================================")

d = np.float16(1)/np.float16(9801)
print("numpy float16   : {:.50f}".format(d))

e = np.float32(1)/np.float32(9801)
print("numpy float32   : {:.50f}".format(e))

f= np.float64(1)/np.float64(9801)
print("numpy float64   : {:.50f}".format(f))

print("========================================================================")

print("float(default)  : {:.50f}".format(1/9801))

> numpy half      : 0.00010204315185546875000000000000000000000000000000
> numpy float     : 0.00010203040506070808982747033510918299725744873285
> numpy double    : 0.00010203040506070808982747033510918299725744873285
> ========================================================================
> numpy float16   : 0.00010204315185546875000000000000000000000000000000
> numpy float32   : 0.00010203040437772870063781738281250000000000000000
> numpy float64   : 0.00010203040506070808982747033510918299725744873285
> ========================================================================
> float(default)  : 0.00010203040506070808982747033510918299725744873285

というふうにデフォルトのfloatがすでに64bitの豪華な型だったのでこれ以上の精度は出せませんでした。

一応 np.longdouble および np.float128 なる型が存在するのですが、これは割り当てられるビット数がハードウェアやOSに依存するということで、paiza.ioやローカルの環境では64bitに抑え込まれてしまいました。

(64bitより大きい型を実現できる環境を作りたかったのですが、断念してしまいました。)

print("float(default)  : {:.50f}".format(1/9801))

print("========================================================================")

h = np.longdouble(1)/np.longdouble(9801)
print("numpy longdouble: {:.50f}".format(h))
i = np.float128(1)/np.float128(9801)
print("numpy float128  : {:.50f}".format(i))

> float(default)  : 0.00010203040506070808982747033510918299725744873285
> ========================================================================
> numpy longdouble: 0.00010203040506070808982747033510918299725744873285
> numpy float128  : 0.00010203040506070808982747033510918299725744873285

numpyの数値型に関する詳細はこちらで。

Data types — NumPy v1.17 Manual

終わりに

結論としては、Decimalを使うのが最強だけど、日常生活でそこまで気にすることはないだろうという感想です。
いつ使うねんこんな知識という感じになってしまいましたが、なんだかんだ色々と知らないことが知れたので良かったです。
ここまでお読みいただきありがとうございました。

*1:ループ、goto文、再帰を用いずに1-100を出力するには？ https://www.quora.com/How-can-I-print-1-to-100-in-C++-without-a-loop-goto-or-recursion

*2:環境の差異とかが面倒なので、すべてpaiza.io上での実行結果を表示しています。

*3: 1/9801 = 0.000102030405060708091011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969799000

*4:2020/04の言語アップデートによりmathモジュールが使えるようになりました。AtCoderさんありがとうございます！