プロジェクトオイラー
http://projecteuler.net/index.php
Q212.
与えられた乱数で直方体を50000個生成する。その占める体積を求めよ。
50000個の直方体同士の重なりを計算すると大変なので、よくやるように小さなセルに区切ってそこに含まれる直方体同士で重なりの計算をする。ここでは1辺の長さ400の立方体に区切る。
例によって、最初に直方体の体積の和を取り、偶数個の直方体の重なりはマイナスし、奇数個の重なりはプラスする。
3個以上の重なりの判定は次のように行う。まず2個の重なりを計算する。そのとき、重なった直方体同士を線で結んでグラフを作る。そして、その部分グラフが3点の完全グラフならその3つの直方体が重なりを持つ可能性があるので、重なりの体積の計算をする。
from itertools import count, product, combinationsdef gen_random():
s = [ 0 ]
for k in range(1, 56):
n = int(((300007 * k * k - 200003) * k + 100003) % 1000000)
yield n
s.append(n)
for k in count(56):
n = (s[(k-24)%56] + s[(k-55)%56]) % 1000000
yield n
s[k%56] = ndef gen_cuboid():
g = gen_random()
while True:
yield (g.next() % 10000, g.next() % 10000, g.next() % 10000,
g.next() % 399 + 1, g.next() % 399 + 1, g.next() % 399 + 1)def intersect_core(x1, dx1, x2, dx2):
if x1 <= x2 and x2 < x1 + dx1:
return (x2, min(x1 + dx1, x2 + dx2) - x2)
elif x2 <= x1 and x1 < x2 + dx2:
return (x1, min(x2 + dx2, x1 + dx1) - x1)
else:
return (0, 0)def intersect(c1, c2):
x, dx = intersect_core(c1[0], c1[3], c2[0], c2[3])
if dx == 0:
return 0
y, dy = intersect_core(c1[1], c1[4], c2[1], c2[4])
if dy == 0:
return 0
z, dz = intersect_core(c1[2], c1[5], c2[2], c2[5])
if dx == 0:
return 0
return (x, y, z, dx, dy, dz)def volume(c):
return c[3] * c[4] * c[5]def is_neighbor(c1, c2, g):
for c in g[c1]:
if c == c2:
return True
return Falsedef is_perfect_sub_graph(g, c, e):
for c1 in e:
for c2 in e:
if c1 == c2:
break
if not is_neighbor(c1, c2, g):
return False
return Truedef intersect_volume(a):
for i in range(len(a) - 1):
if i == 0:
c = intersect(a[0], a[1])
else:
c = intersect(c, a[i+1])
if c == 0:
return 0
return volume(c) * (-1) ** (len(a) - 1)def add_graph_core(g, c1, c2):
if c1 in g:
g[c1].append(c2)
else:
g[c1] = [ c2 ]def add_graph(g, c1, c2):
add_graph_core(g, c1, c2)
add_graph_core(g, c2, c1)class space:
def __init__(self, s, l):
self.span = s
self.cell_length = l
self.num = s / l
self.cells = [ [ ] for i in range(self.num ** 3) ]
self.cuboids = [ ]
def add(self, c):
id = len(self.cuboids)
self.cuboids.append(c)
x0 = c[0] / self.cell_length
y0 = c[1] / self.cell_length
z0 = c[2] / self.cell_length
x1 = (c[0] + c[3] - 1) / self.cell_length
y1 = (c[1] + c[4] - 1) / self.cell_length
z1 = (c[2] + c[5] - 1) / self.cell_length
for x, y, z in product(range(x0, x1 + 1),
range(y0, y1 + 1), range(z0, z1 + 1)):
cell_id = (z * self.num + y) * self.num + x
self.cells[cell_id].append(id)
def volume(self):
total_vol = sum(map(volume, self.cuboids))
set_intersection = set()
for cell in self.cells:
if len(cell) > 1:
for s, v in self.gen_intersection(cell):
if s not in set_intersection:
total_vol += v
set_intersection.add(s)
return total_vol
def gen_intersection(self, cell):
g = { }
for i in range(len(cell)):
i1 = cell[i]
c1 = self.cuboids[i1]
for j in range(i + 1, len(cell)):
i2 = cell[j]
c2 = self.cuboids[i2]
c3 = intersect(c1, c2)
if c3:
add_graph(g, i1, i2)
yield ( (i1, i2), -volume(c3) )
for n in count(3):
exists = False
for s, v in self.gen_multi_intersection(g, n):
exists = True
yield (s, v)
if not exists:
break
def gen_multi_intersection(self, g, n):
for i, a in g.iteritems():
if len(a) < n - 1:
continue
for e in combinations(a, n - 1):
if e[0] <= i:
continue
if is_perfect_sub_graph(g, i, e):
t = (i,) + e
u = tuple(map(lambda k: self.cuboids[k], t))
yield (t, intersect_volume(u))N = 50000
S = 10400
L = 400 # edge length of cell
s = space(S, L)
g = gen_cuboid()
for n in range(N):
s.add(g.next())
print s.volume()
