次の関数について考えます。
これをそのまま書くと、kはnまでで十分だから、
def H(n): return sum(n / k for k in xrange(1, n + 1))
ですね。O(n)の計算量のはずです（nが小さければ）。
さて、この計算を速くできないでしょうか。まず、n = 10で具体的な計算をしてみましょう。
H(n) = 10/1 + 10/2 + 10/3 + 10/4 + 10/5 + 10/6 + 10/7 + 10/8 + 10/9 + 10/10
= 10 + 5 + 3 + 2 + 2 + 1 + 1 + 1 + 1 + 1
です。最後のほう、1が続きますよね。n / 2個続くはずです。ここを一々割り算しているのはムダですよね。1になるkはn/2 + 1 〜 nまでです。lになるkはn/(l+1) + 1 〜 n/lで、このように商がkになる個数はO(1)で求められます。
kが√nまでは素直に求めて、それ以降はこの方法で求めましょう。最大の商はn/√n = √nだから、どちらも計算量はO(√n)で、合わせてもO(√n)となります。
from math import sqrt def H(n): m = int(sqrt(n)) s = sum(n / k for k in xrange(1, m + 1)) start = n for k in xrange(1, n / m): end = start start = n / (k + 1) s += k * (end - start) return s
このように領域を分けると計算法を変えると計算量が劇的に抑えられることがよくあります。特に上の例のように√nで分けるとうまくいくことがProject Eulerではときどきあります。
余談ですが、上のH(n)はよく考えると、
{ x, y | xy <= n, x > 0, y > 0 }
という集合の要素数であることがわかります。x > yだけ考えて2倍して、別途x = yを計算すると、上の方法より2倍以上速くなりました。
from math import sqrt def H(n): m = int(sqrt(n)) return sum(n / k - k for k in xrange(1, m + 1)) * 2 + m
