Project Euler 210

プロジェクトオイラー
http://projecteuler.net/index.php?section=problems&id=210

これにもやはり数値計算上のトラップが。

この問題は、おおざっぱに言うと、円の内側にある格子点の個数を数える。

例えば、xy平面を考えて、中心が原点で半径rが10だとする。そのとき第1象限におけるy=7上の格子点の数は、√(100-49) = 7.14…だから、7つ。y=8なら、√(100-64) = 6だから、円周上の格子点は数えないので、5つ。これをまとめると、数学的には、

[√(r² - y² - 0.5)]

個とすればよい。

しかし、IEEE754の倍精度は15.6〜16桁程度しか精度がないので、整数が大きくなってくると0.5を引く手が使えなくなってくる。0.5を引いても値が変わらないのだ。例えば、半径が10⁹で、y=8×10⁸上の格子点の個数は、

from math import sqrt
def num_of_lattice_points(r, y):

	return int(sqrt(r * r - y * y - 0.5))
R = 10 ** 9

Y = 8 * 10 ** 8

print num_of_lattice_points(R, Y)	# 600000000

となって、正しい答えより1個多くなってしまう。この場合、もう一度整数計算に戻って確認する必要がある。あるいは、sqrtの計算も整数計算にするか。

この問題は、このトラップが使えるギリギリの半径が与えられていた。