Project Euler 126

http://projecteuler.net/index.php?section=problems&id=126

まず、a × b × cでm層目の立方体の数は、詳しくは書きませんが、サイドとそれ以外を分けて考えると易しいです。

P ≡ a + b + c
Q ≡ b c + c a + a b

とすると、

A_m = 4(m-1)(P+m-2) + 2Q

となります。
さて、この問題はA_mを小さい順に並べにくいので難しいです。ある数以下というのは比較的易しいです。a, b, cが決まっていれば1層目が最も小さいので、2(b c + c a + a b)がある数以下の組み合わせについて調べればよいです。そのため、そのある数をだんだん大きくしていけばいつか答えに辿り着くでしょう。
このだんだん大きくするというのはどうするのがよいでしょう。これはなかなか難しいと思われます。この問題はO(n²)程度のようです。2倍ずつしていくことにしましょう。例えば1000から始めて2000、4000とします。そうすると、例えば16001が答えだと32000の計算をしなければなりません。すなわち最も適切な数から始めたときと比べて4倍くらいの計算量になります。もっと倍率を小さくすればよいですが、そうすると重複の計算が増えてしまいます。O(n^e)として増やしていく倍率をaとすると、一発で答えを出すときに比べてa^e / (a^e - 1)倍の計算量になります。さきほどの最悪の場合の計算量をかけると、a^2e / (a^e - 1)となります。最悪の場合を最小にすることを考えると、

a^2e / (a^e - 1) = 1 / (a^-e - a^-2e)

この分母が最大になるのはa^-e = 1 / 2だから、

a = 2^1/e

この問題は、1.4倍としてみました。

from itertools import *

def iterate(f, x):
    while True:
        yield x
        x = f(x)

def gen_PQ(max_Q):
    for a in xrange(1, (max_Q - 1) / 2 + 1):
        for b in xrange(1, min(a, (max_Q - a) / (a + 1)) + 1):
            for c in xrange(1, min(b, (max_Q - a * b) / (a + b)) + 1):
                yield a + b + c, b * c + c * a + a * b

def calc_C(max_n):
    C = [ 0 ] * (max_n + 1)
    for P, Q in gen_PQ(max_n / 2):
        for n in takewhile(lambda n: n <= max_n,
                (4 * (m - 1) * (P + m - 2) + 2 * Q for m in count(1))):
            C[n] += 1
    return C

def find(C, N):
    try:
        return C.index(N)
    except:
        return -1

N = 1000
print next(n for n in (find(calc_C(max_n), N) for max_n in
                    iterate(lambda x: int(x * 1.4), 1000)) if n != -1)