プロジェクトオイラー
http://projecteuler.net/index.php
Q202.
正三角形状に鏡を配置する。頂点をA,B,Cとする。その頂点には無限小の隙間があり、Cからレーザーが入って反射する。12017639147回レーザーが反射して再びCから出る入射法は何通りあるか。
よくあるように、鏡に当たっても反射しないで仮想的に通過すると考えるとよい。そうすると、鏡は正三角形が無限に積み重なるように配置され、それを通過することになる。レーザーがある頂点にたどり着いたとすると、その座標は、
sCA + tCB
と表される。sとtは整数である。s-tが3の倍数なら、これはCとなる。
ここにたどり着くまでに、2(s + t) - 3回鏡に当たっている。だから、N回反射してCに戻ってくるとすると、n = s + t = (N + 3) / 2で表される直線上の頂点にたどり着くことになる。そのような頂点はn - 1個あり、その中でs - tが3の倍数のものを選べばよい。ただし、sとtが互いに素でなければ、その前にどこかの頂点から出ていることになる。nを素因数分解して、でてきた素数の倍数にsがなったときは排除する。排除の仕方は、例によって和集合の公式を使う。
from math import sqrtdef is_prime(n):
for p in primes:
if p * p > n:
return True
elif n % p == 0:
return False
return Truedef make_primes(n):
m = int(sqrt(n + 0.5))
for p in xrange(3, m + 1, 2):
if is_prime(p):
primes.append(p)
a = [ 0 ] * ((n + 15) / 16)
for p in primes:
if p * p > n:
break
if p == 2:
continue
for k in xrange(3 * p, n, 2 * p):
a[k>>4] |= 1 << (k & 15)
for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):
if (a[k>>4] & (1 << (k & 15))) == 0:
primes.append(k)def factorize(n, q = 2):
if n == 1:
return [ ]
elif n < q * q:
return [ n ]
for p in primes:
if p < q:
continue
if n % p == 0:
return [ p ] + factorize(n / p, p)
return [ n ]def unique(a):
i = 0
while i + 1 < len(a):
if a[i] == a[i+1]:
a.pop(i)
i += 1def gen_free_square(n):
a = factorize(n)
unique(a)
for e in gen_free_square_core(a):
yield edef gen_free_square_core(a, i = 0):
yield 1
if i < len(a):
for k in range(i, len(a)):
p = a[k]
for e in gen_free_square_core(a, k + 1):
yield -p * edef num_solution(n, e):
positive = e > 0
e = abs(e)
if n % 3 == 1:
# (3k + 2, n - 3k - 2)
x0 = (3 - e % 3) * e
k0 = (x0 - 2) / 3
kmax = (n - 2) / 3
m = (kmax - k0) / e + 1
elif n % 3 == 2:
# (3k + 1, n - 3k - 1)
x0 = e % 3 * e
k0 = (x0 - 1) / 3
kmax = (n - 2) / 3
m = (kmax - k0) / e + 1
else:
m = 0
if positive:
return m
else:
return -mN = 12017639147
n = (N + 3) / 2
primes = [ 2 ]
make_primes(int( (n + 0.5) ** 0.5))s = 0
for e in gen_free_square(n):
s += num_solution(n, e)
print s