Pythonでは整数は大きくなると勝手に多倍長整数として表すので誤差は出ませんが、浮動小数点数を介すると途端に誤差が発生します。Project Eulerで特に問題になるのは、整数の平方根の計算です。すなわち、
⌊√n⌋
の計算です。この計算は整数nを浮動小数点数に変換して平方根の計算をします。浮動小数点数は15桁くらい精度があるため、それくらいまでの大きさの整数なら、
int(sqrt(n))
で正確な計算をしてくれます。しかし、それより大きい整数だとそうもいきません。例えば、1234567892098765184は19桁の整数ですが、
n = 1234567892098765184 m = int(sqrt(n)) print m ** 2 - n # 265
とmは√nより大きくなってしまいます。これは、整数1234567892098765184は浮動小数点数化すると精度が足りず、1234567892098765312.0と大きくなってしまうからです。
⌊√n⌋
を精度よく計算するには、まず精度を上げることです。decimalモジュールを使うと精度を自由に変えることができます。
from decimal import * getcontext().prec = 20 # デフォルト28桁のところを20桁にする m = int(Decimal(n).sqrt()) print m ** 2 - n # -2222221948
ただ、精度が上がる分遅くなります。
整数の範囲で平方根を計算することもできます。
def int_sqrt(n): def f(prev): while True: m = (prev + n / prev) / 2 if m >= prev: return prev prev = m return f(n)
m <- (n + n / m) / 2
を繰り返して収束させる方法です。ただ、Pythonでこれは遅いです。といっても、上の精度を上げる方法より少し速いようです。
この方法は、初期値によって速度が変わります。例えば√30を考えましょう。30から始めると、
30 -> (30 + 30/30) / 2 = 15 -> (15 + 30/15) / 2 = 8 -> (8 + 30/8) / 2 = 5 -> (5 + 30/5) / 2 = 5
で5となりますが、6から始めると、
6 -> (6 + 30/6) / 2 = 5 -> (5 + 30/6) / 2 = 5
で収束します。このように当然ながら√nに近い初期値なら速くなります。初期値は真値以上でなければなりません。だから、int(sqrt(n))より少し大きいとよいです。int(sqrt(n) * (1 + 1e-10))くらいでよいでしょう。
def int_sqrt2(n): def f(prev): while True: m = (prev + n / prev) / 2 if m >= prev: return prev prev = m return f(int(sqrt(n) * (1 + 1e-10)))
こうすると速くなります。10万回回すと、Pythonで順番に1.3s, 2.5s, 0.2s、PyPyで0.8s, 0.9s, 0.1sでした。
