Project Euler 95

Problem 95
真約数はその数自信を除いた約数全てである。例えば、28の真約数は1,2,4,7,14となる。その和は28となるので、これを完全数と呼ぶ。
面白いことに220の真約数の和は284で284の真分数の和は220であり、これは2つの数から成るチェーンを成す。これより220と284は友愛数と呼ぶ。
恐らくさらに長いチェーンはあまり知られていない。例えば、12496から始めると、5つの数から成るチェーンを形成する。
12496 → 14288 → 15472 → 14536 → 14264 (→ 12496 → ...)
このチェーンは開始点に戻るので、友愛鎖と呼ばれる。
どの要素も100万を超えない最も長い友愛鎖の最小数を求めよ。
http://projecteuler.net/index.php?section=problems&id=95

あらかじめ真約数の和をエラトステネスのふるい的に求めておきます。そして、小さいほうから順番に1に辿り着くかループに辿り着くかループになるのかを調べます。1に辿り着くかループに辿り着いたら0を記録し、自分に戻ってきたらループの長さを記録します。最後に最もループが長くて最小の数を得ます。

from itertools import ifilter, izip
from math import sqrt

def prime_pow_sum(n, p):
    s = 1
    while n % p == 0:
        s = s * p + 1
        n /= p
    return s, n

def calc_sum_divisors(max_n):
    a = range(max_n + 1)
    b = [ 1 ] * (max_n + 1)
    for p in ifilter(lambda p: a[p] == p, xrange(2, int(sqrt(max_n)) + 1)):
        for n in xrange(p, max_n + 1, p):
            s, m = prime_pow_sum(a[n], p)
            a[n] = m
            b[n] *= s
    
    for n in xrange(2, max_n + 1):
        if a[n] != 1:
            b[n] *= (a[n] + 1)
        b[n] -= n
    
    return b

def calc_loop_length(n):
    if b[n] >= 0:
        return
    
    c = [ n ]
    s = set(c)
    while True:
        p = c[-1]
        q = a[p]
        
        if q > N:               # OB
            for m in c:
                b[m] = 0
            return
        elif q == n:            # loop
            for m in c:
                b[m] = len(c)
            return
        elif b[q] >= 0:         # to 1 or branch
            for m in c:
                b[m] = 0
            return
        elif q in s:            # branch of new loop
            k = c.index(q)
            for m in (c[i] for i in xrange(k)):
                b[m] = 0
            for m in (c[i] for i in xrange(k, len(c))):
                b[m] = len(c) - k
            return
        
        c.append(q)
        s.add(q)


N = 10 ** 6
a = calc_sum_divisors(N)
b = [ -1 ] * (N + 1)
b[1] = 0
for n in xrange(2, N + 1):
    calc_loop_length(n)
def longer(x, y): return x if x[1] >= y[1] else y
print reduce(longer, izip(xrange(N + 1), b))