Project Euler #23 - Non-abundant sums

Non-abundant sums: 두 개의 abundant numbers(자기 자신을 뺀 약수의 합이 자기 자신 보다 커지는 숫자)를 더해서 만들 수 없는 숫자들의 합은?

문제에서 28123보다 큰 수는 모두 두 개의 a.n 합으로 표현할 수 있다고 했으니 그보다 작은 수만 체크하면 된다. 두 개의 a.n 로 만들 수 "없는" 것보다는 만들 수 "있는" 쪽을 찾는 게 쉬워 보인다.

sympy를 쓰면 어렵지 않다.

from sympy import divisors

abun = []
for i in range(12,28123):
    if sum(divisors(i)[:-1]) > i: abun.append(i)

allabuncomb = [True] * 28123
for i,a in enumerate(abun):
    for b in abun[i:]:
        if a+b >= 28123: break
        allabuncomb[a+b] = False

print sum([n for n in range(28123) if allabuncomb[n]])

역시나 다른 문제들처럼 소수-소인수분해-약수합 부분을 직접 구현하면 약간 빨라진다. 하지만 abundant number 를 찾는 부분보다, 그 숫자들로 조합을 만드는 부분이 오래 걸리는 거라 큰 개선은 없다.

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 isAbun(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 > n)

pn = rwh_primes1(28123)
abun = []
for n in range(12,28123):
    if isAbun(n,pn): abun.append(n)

allabuncomb = [True] * 28123
for i,a in enumerate(abun):
    for b in abun[i:]:
        if a+b >= 28123: break
        allabuncomb[a+b] = False

print sum([n for n in range(28123) if allabuncomb[n]])