Project Euler #21 - amicable numbers

10000이하의 친화수를 찾는 문제

sympy의 divisors함수를 쓰면 쉽다.

from sympy import divisors

d_nDic = {}
for i in range(3,10000):
    d_nDic[i] = sum(divisors(i))-i

ami = []
for i in d_nDic:
    if d_nDic[i] > 3 and d_nDic[i] < 10000:
        if i == d_nDic[d_nDic[i]] and i != d_nDic[i]:
            ami.append(i)

print sum(ami)

10000까지 반복해서 소인수분해를 해야 하는데, 좀 더 빠르게 하려면 소수를 직접 찾고, 소인수 분해를 직접 한 다음에, 약수의 합을 계산하면 된다.
마지막 스텝은 중학교 때 배운 공식 24 = 2^3 * 3 => 24의 약수의 합=(1+2+4+8)*(1+3)=1+2+4+8+3+6+12+24 를 쓰면 된다.
코드는 길어지지만 실행시간은 1/4로 줄어든다.

def rwh_primes1(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

def factor(n,pn):
    rs = {}
    for p in pn:
        if n % p == 0:
            cnt = 0
            while n % p == 0:
                n /= p
                cnt += 1
            rs[p] = cnt
            if n == 1: return rs
        if p*p > n:
            rs[n] = 1
            return rs

def d_n(n,pn):
    fDic = factor(n,pn)
    dn = 1
    for k in fDic:
        dn *= sum([k**i for i in range(fDic[k]+1)])
    return dn - n

pn = rwh_primes1(10000)
d_nDic = {}
for i in range(3,10000):
    d_nDic[i] = d_n(i,pn)

ami = []
for i in d_nDic:
    if d_nDic[i] > 3 and d_nDic[i] < 10000:
        if i == d_nDic[d_nDic[i]] and i != d_nDic[i]:
            ami.append(i)

print sum(ami)