Problem 86
蜘蛛Sが直方体で大きさが6×5×3の部屋の一つの角に座っており、蝿Fが反対側の角に座っている。部屋の表面を進んでSからFの最短の「直線」距離は10（中略）。
どの与えられた直方体にも「最短」経路候補が3つまであり、その最短ルートはいつも整数であるとは限らない。
大きさが整数でM×M×Mまでの直方体の部屋すべてを考えると、M=100のとき最短距離が整数になる直方体は2060ある。そしてこれは解の個数が2000を最初に超える最小のMである。M=99なら解の個数は1975である。
解の個数が100万を最初に超える最小のMを求めよ。
http://projecteuler.net/index.php?section=problems&id=86
直方体の辺の長さをp ≤ q ≤ rとすると、最短距離は
√((p + q)2 + r2)
となります。要するに、p + qとrとこれで直角三角形を成します。
いつものようにピタゴラス数の生成式を使います。
a = l(m2 - n2) b = 2lmn c = l(m2 + n2)
例えば大きさ6までとして、(6, 8, 10)という直角三角形が生成されたとき、8を分割して
(2, 6, 6), (3, 5, 6), (4, 4, 6)
という直方体が生成されます。大きさが8までなら6,8のどちらかを分割して、
(1, 5, 8), (2, 4, 8), (3, 3, 8), (2, 6, 6), (3, 5, 6), (4, 4, 6)
という直方体が距離10になります。このようにして直方体の数を数えて、Mを一つずつ大きくして直方体の数が100万を超えたところでやめます（Code1）。
しかしこれは遅いです。Mを探索するのに1ずつ増すのではなく、二分探索を使いましょう。その前に1から始めて倍々して100万を超えるMを求め、そのMと前のMの間で二分探索します。これでかなり速くなります（Code2）。
より速い方法があります。それは少なくとも一つの辺がちょうどMになる直方体を求めるものです。それには直角三角形の斜辺以外のどちらかがMになる必要があります。すなわち、ピタゴラス数の生成式のl,m,nを使ってM = l(m2 - n2)またはM = 2lmnにならなければなりません。そこで、Mを素因数分解して約数を求めて個数を数えていきます（Code3）。
素因数分解をエラトステネスのふるい的にしてもみましたが、あまり速くなりませんでした。
# Code1 from itertools import imap, ifilter, count from fractions import gcd from math import sqrt def head(a): for e in a: return e def S(M): def f(m): counter = 0 for n in ifilter(lambda n: gcd(m, n) == 1, xrange(2 if m & 1 else 1, m, 2)): a = m * m - n * n b = 2 * m * n if a > b: a, b = b, a for l in xrange(1, M / a + 1): c, d = a * l, b * l if d <= M: counter += c / 2 if c * 2 > d: counter += c - (d - 1) / 2 elif c <= M and d < c * 2: counter += (c * 2 - d + 2) / 2 return counter return sum(imap(f, xrange(2, (M + 1) / 2 + 1))) N = 10 ** 6 a = head(ifilter(lambda x: x[1] > N, ((M, S(M)) for M in count(1)))) print a[0]
# Code2 from itertools import imap, ifilter from fractions import gcd from math import sqrt # SはCode1と同じ def bin_search(begin, end): if end - begin == 1: return begin else: mid = (begin + end - 1) / 2 n = S(mid) if n <= N: return bin_search(mid + 1, end) else: return bin_search(begin, mid + 1) N = 10 ** 6 M = 1 while True: n = S(M) if n > N: break M *= 2 print bin_search(M / 2 + 1, M + 1)
# Code3 from itertools import imap, ifilter, count, takewhile from fractions import gcd from math import sqrt def head(a): for e in a: return e def is_prime(n): return all(n % p != 0 for p in takewhile(lambda p: p * p <= n, primes)) def gen_primes(): for k in count(): if k == len(primes): primes.append(head(ifilter(is_prime, count(primes[-1] + 1)))) yield primes[k] def div_pow(n, p): e = 0 while n % p == 0: n /= p e += 1 return n, e def factorize(n): if n == 1: return [ ] for p in takewhile(lambda p: p * p <= n, gen_primes()): n, e = div_pow(n, p) if e > 0: return [(p, e)] + factorize(n) return [(n, 1)] def gen_divisors(f, k = 0): if k == len(f): yield 1, 1 else: p, e = f[k] n = p ** e for d1, d2 in gen_divisors(f, k + 1): if p == 2: if e >= 2: yield d1 * n, d2 elif e == 0: yield d1, d2 elif e > 0: yield d1 * n, d2 yield d1, d2 * n else: yield d1, d2 def gen_divisors2(f, k = 0): if k == len(f): yield 1, [] else: p, e0 = f[k] for d1, f2 in gen_divisors2(f, k + 1): for e in range(e0 + 1): yield d1 * p ** e, [(p, e0 - e)] + f2 def gen_divs(f): for d1, f2 in gen_divisors2(f): for d2, d3 in gen_divisors(f2): yield d1, d2, d3 def count_cuboid(M, b): if b < M: return b / 2 elif b < M * 2: return M - (b - 1) / 2 else: return 0 def T(M): counter = 0 f = factorize(M) for l, d1, d2 in gen_divs(f): if (d1 + d2) % 2 == 0: if d1 <= d2: continue m = (d1 + d2) / 2 n = m - d2 b = 2 * l * m * n else: m = d1 / 2 n = d2 if m < n: m, n = n, m b = l * (m * m - n * n) counter += count_cuboid(M, b) return counter def gen_S(): s = 0 for M in count(1): s += T(M) yield M, s N = 10 ** 6 primes = [ 2 ] a = head(ifilter(lambda x: x[1] > N, gen_S())) print a[0]
