プロジェクトオイラー
http://projecteuler.net/index.php
Q199.
図のように3つの円に接するように円を描いていく。これを10回繰り返したとき、円で塗りつぶされていない領域の割合はいくらか？
互いに接する3つの円C1,C2,C3の半径をそれぞれr1,r2,r3、その中心をO1,O2,O3、aをO1からO2へのベクトル、bをO1からO2へのベクトル、3つの円の間に描いた円の半径をrとすると、O1からその円の中心へのベクトルはs,tを使って、
sa + tb
と表せる。
そうすると、次の式が成り立つ。
sa + tb = r + r1 (s-1)a + tb = r + r2 sa + (t-1)b = r + r3
この3式を各々自乗して、展開する。
これを解けばいいのだが、いくらやっても解けなかった。仕方がなく、3変数のニュートン法を使って解いた。
これが解ければあとは再帰で簡単。
from math import sqrtdef calc_R(r1, r2, r3):
eps = 1e-8
a2 = (r1 + r2) ** 2
b2 = (r1 + r3) ** 2
ab = r1 * (r1 + r2 + r3) - r2 * r3
def dist(s_, t_, r_):
return a2 * s_ ** 2 + 2 * ab * s_ * t_ + b2 * t_ ** 2 - r_ ** 2
def f(s, t, r):
return dist(s, t, r + r1)
def f_s(s, t, r):
return 2 * (a2 * s + ab * t)
def f_t(s, t, r):
return 2 * (ab * s + b2 * t)
def f_r(r):
return -2 * (r + r1)
def g(s, t, r):
return dist(s - 1, t, r + r2)
def g_s(s, t, r):
return 2 * (a2 * (s - 1) + ab * t)
def g_t(s, t, r):
return 2 * (ab * (s - 1) + b2 * t)
def g_r(r):
return -2 * (r + r2)
def h(s, t, r):
return dist(s, t - 1, r + r3)
def h_s(s, t, r):
return 2 * (a2 * s + ab * (t - 1))
def h_t(s, t, r):
return 2 * (ab * s + b2 * (t - 1))
def h_r(r):
return -2 * (r + r3)
def is_diff_large(r1, r2, r3):
r_max = max(r1, r2, r3)
if r_max == r1:
if r2 < 0.5 * r1 and r3 < 0.5 * r1:
return 1
else:
return 0
elif r_max == r2:
if r3 < 0.5 * r2 and r1 < 0.5 * r2:
return 2
else:
return 0
else:
if r1 < 0.5 * r3 and r2 < 0.5 * r3:
return 3
else:
return 0
def init_st():
if r1 > 0 and r2 > 0 and r3 > 0:
mode = is_diff_large(r1, r2, r3)
if mode > 0:
if mode == 1:
r12, r22, r32 = r1, r2, r3
elif mode == 2:
r12, r22, r32 = r2, r3, r1
else:
r12, r22, r32 = r3, r1, r2
r2r = sqrt(r22 * r12)
r3r = sqrt(r32 * r12)
m = (1 - sqrt(r22 * r22 + r32 * r32)
/ (r12 * 2)) / (r2r + r3r)
s = r3r * m
t = r2r * m
if mode == 2:
s, t = 1 - s - t, s
elif mode == 3:
s, t = t, 1 - s - t
else:
r23 = (r2 + r3) ** 2
r31 = (r3 + r1) ** 2
r12 = (r1 + r2) ** 2
r_all = r23 + r31 + r12
s = r31 / r_all
t = r12 / r_all
else:
r2r = sqrt(-r2 * r1)
r3r = sqrt(-r3 * r1)
m = (1 - sqrt(r2 * r2 + r3 * r3)
/ (r1 * 2)) / (r2r + r3r)
s = r3r * m
t = r2r * m
return (s, t)
def op(v1, v2):
return (v1[1] * v2[2] - v2[1] * v1[2],
v1[2] * v2[0] - v2[2] * v1[0],
v1[0] * v2[1] - v2[0] * v1[1])
def ip(v1, v2):
return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
s, t = init_st()
r = 0.0
while True:
va = (f_s(s, t, r), f_t(s, t, r), f_r(r))
vb = (g_s(s, t, r), g_t(s, t, r), g_r(r))
vc = (h_s(s, t, r), h_t(s, t, r), h_r(r))
vf = (-f(s, t, r), -g(s, t, r), -h(s, t, r))
vbc = op(vb, vc)
vca = op(vc, va)
vab = op(va, vb)
abc = ip(va, vbc)
ds = (vbc[0] * vf[0] + vca[0] * vf[1] + vab[0] * vf[2]) / abc
dt = (vbc[1] * vf[0] + vca[1] * vf[1] + vab[1] * vf[2]) / abc
dr = (vbc[2] * vf[0] + vca[2] * vf[1] + vab[2] * vf[2]) / abc
if abs(ds) + abs(dt) + abs(dr) < eps:
return r + dr
s += ds
t += dt
r += drdef cover_three_circles_gap(r1, r2, r3, n):
if n == 0:
return 0
r = calc_R(r1, r2, r3)
area1 = cover_three_circles_gap(r, r2, r3, n - 1)
area2 = cover_three_circles_gap(r1, r, r3, n - 1)
area3 = cover_three_circles_gap(r1, r2, r, n - 1)
return r * r + area1 + area2 + area3N = 10
orig_area = 1
r = 3 / (3 + 2 * sqrt(3))
cover_area = 3 * r ** 2 + 3 * cover_three_circles_gap(-1, r, r, N)
cover_area += cover_three_circles_gap(r, r, r, N)
print "%.8f" % (1 - cover_area / orig_area, )
