プロジェクトオイラー
http://projecteuler.net/index.php
Q184.
中心が原点、半径rの円の内側（周上は含まない）の格子点からなる三角形で原点を内部に含むものの個数をIrとする。I105を求めよ。
対称性を考えると、三角形の点の位置は次の5種類に分類される。
i) 第3象限、第3象限、第1象限
ただし、2番目の点は1番目の点と原点を結ぶ直線より右下にある
ii) 第3象限、第4象限、第1象限
iii) (x, 0)(x < 0)、第3象限、第1象限
iv) (x, 0)(x < 0)、第4象限、第1象限
v) (x, 0)(x < 0)、(0, y)(y < 0)、第1象限
このうち、iiiのみ反転と回転で1つの三角形から8つが生成される。それ以外は回転のみの4つである。
計算するときに、その点と原点を結ぶ直線より左上の第1象限の点の個数を返す関数fを作っておくと便利である。例えばiiなら、2番目の点は自由に選べて、1番目の点を(-x, -y)とすると、3番目の点はf(x, y, r)個選べる。iについては、3番目の点を(x, y)に選ぶと、1番目の点はf(y, x, r)個、2番目の点はf(x, y, r)個それぞれ独立に選べる。
この関数は度々使われるので、あらかじめ配列にしておくと速くなる。そうすると、この配列を作る時間が全体のほとんどを占めるようになる。そのため、この配列を速くする。それには、第1象限の点で、xとyが互いに素のものを選び、同じ比になる組合せの個数を配列にしておく。そして、それをその点と原点を結ぶ直線とx軸が成す角の大きいほうから順に並べ替えておく。そうすると、fの配列の計算はただの足し算になり、計算は一瞬で終わる。
from math import sqrt, ceil
import fractionsdef gen_pt(r):
for y in xrange(1, r):
xmax = int(ceil(sqrt(r * r - y * y - 1e-9)))
for x in xrange(1, xmax):
yield (x, y)def calc_pt(r):
a = [ ]
for x, y in gen_pt(r):
if fractions.gcd(x, y) == 1:
dmax = int(ceil(sqrt(float(r * r)
/ (x * x + y * y) - 1e-9))) - 1
a.append( (x, y, dmax) )
a.sort(cmp = lambda pt1, pt2: pt1[0] * pt2[1] - pt1[1] * pt2[0])
return adef calc_f_core(x1, y1, r):
n = 0
for x2, y2 in gen_pt(r):
if x1 * y2 - x2 * y1 > 0:
n += 1
return ndef calc_f(r):
f = [ [ 0 for i in xrange(r) ] for i in xrange(r) ]
n = 0
for x, y, dmax in a:
for d in xrange(1, dmax + 1):
f[d*y][d*x] = n
n += dmax
n = calc_f_core(1, 0, r)
for x in xrange(1, r):
f[0][x] = n
for y in xrange(1, r):
f[y][0] = 0
return fdef g(x, y, r):
d = fractions.gcd(x, y)
x /= d
y /= d
return int(ceil(sqrt(r * r / (x * x + y * y) - 1e-9)))def calc(r):
counter = 0
# i)
for x1, y1 in gen_pt(r):
counter += f[y1][x1] * f[x1][y1]
# ii)
for x1, y1 in gen_pt(r):
counter += f[0][1] * f[y1][x1]
# iii)
for x1, y1 in gen_pt(r):
counter += (r - 1) * f[x1][y1] * 2
counter += (r - 1) * f[0][1] ** 2 # iv)
counter += (r - 1) ** 2 * f[0][1] # v)
counter *= 4
return counterr = 105
a = calc_pt(r)
f = calc_f(r)
print calc(r)
